---
title: "課題3：nested logitの市場レベル需要推定と価格内生性"
subtitle: "計量経済学II"
author: "Kei Ikegami"
lang: ja
format:
  html:
    toc: true
    toc-depth: 3
    toc-location: left
    number-sections: false
    embed-resources: true
    code-copy: true
    theme: cosmo
    fig-width: 8
    fig-height: 4.8
execute:
  warning: false
  message: false
---

```{r}
#| include: false
library(ggplot2)
library(dplyr)
library(fixest)

set.seed(20263)

theme_lecture <- function() {
  fam <- if (Sys.info()[["sysname"]] == "Darwin") "Hiragino Sans" else ""
  theme_minimal(base_size = 13, base_family = fam)
}
theme_set(theme_lecture())
knitr::opts_chunk$set(fig.align = "center", dpi = 144)
options(scipen = 8)
```

::: {.callout-important}
## この課題のゴール

第3回の講義では、缶コーヒー市場（まったり系 vs すっきり系）を例にnested logitとBerry反転、価格内生性、2SLSを学んだ。この課題では、舞台を**炭酸飲料市場**（「甘め系」2財、「無糖・すっきり系」2財）に変え、自分の手でデータを生成し、推定・推論まで一通り実行する。

具体的には次の3つを確認する。

1. plain logitのBerry反転をOLSで推定すると価格係数$\alpha$が過小評価されること。
2. nested logit版のBerry反転をOLSで推定すると、$\alpha$だけでなくネスティング・パラメータ$\rho$も歪むこと。
3. コストシフターと同一ネスト内の競合属性を使った2SLSで、$\alpha$と$\rho$の両方を回収できること。

配点は全体で100点。手を動かせば2〜4時間程度で終わる分量を想定している。ヒントのcalloutを活用してほしい。
:::

## 設定：炭酸飲料市場

大手飲料メーカーの新商品企画チームにいるとしよう。あなたの会社は、全国の$T=100$市場（コンビニチェーンの店舗群、あるいは地域市場と考えてもよい）で、4種類の炭酸飲料を展開している。

- 財1・財2：**「甘め系」**ネスト（フレーバー強め、糖分多め）
- 財3・財4：**「無糖・すっきり系」**ネスト（微炭酸・無糖）
- 財0（outside option）：炭酸飲料を買わない（他のカテゴリの飲み物を買う、あるいは何も買わない）

消費者の誤差項は、講義で扱ったcanonical nested-logitのGEV同時分布に従うとする。この課題では、独立なGumbelを足して誤差を疑似生成せず、講義のwithin/betweenシェア写像をそのままDGPとして使う。平均効用は

$$
\delta_{jt} = x_{jt}'\beta - \alpha p_{jt} + \xi_{jt}
$$

とする。真値は

$$
\alpha = 2.0, \quad \beta = 1.5, \quad \rho = 0.5
$$

とする。

ここで $\lambda=1-\rho$ がdissimilarity parameterである。$0<\lambda\leq1$（$0\leq\rho<1$）が通常のランダム効用モデルと整合的な範囲で、標準的な正規化では同一ネスト内の誤差相関は $1-\lambda^2$。したがって、$\rho$は相関そのものではなく、相関を単調に動かすネスティング・パラメータである。

## Part 1：シミュレーション設計（40点）

### Q1（10点）：外生変数とコストシフターの生成

$T = 100$ 市場、$J = 4$ 財について、次の変数を生成せよ。

- 商品属性 $x_{jt}$：平均3.6、標準偏差1.5の正規分布
- コストシフター $w_{jt}$：平均4.0、標準偏差1.0の正規分布
- 需要ショック $\xi_{jt}$：平均0、標準偏差1.0の正規分布

いずれも市場×財の$100 \times 4$行列として生成すること。`set.seed()`はsetupチャンクで既に固定されているので、以降このチャンクの実行順序を変えないこと。

```{r}
#| label: q1-skeleton
#| eval: false
T_mkt <- 100
J_mkt <- 4

x_mkt <- matrix(rnorm(___, mean = ___, sd = ___), nrow = ___, ncol = ___)
w_mkt <- matrix(rnorm(___, mean = ___, sd = ___), nrow = ___, ncol = ___)
xi_mkt <- matrix(rnorm(___, mean = ___, sd = ___), nrow = ___, ncol = ___)
```

::: {.callout-tip}
## ヒント

`rnorm(n, mean, sd)`は長さ`n`のベクトルを返す。`matrix(rnorm(T*J, ...), nrow = T, ncol = J)`のように書けば、$T \times J$行列に整形できる。第3回講義のRシミュ3と基本的に同じ書き方でよい。
:::

### Q2（15点）：価格の生成（内生性を組み込む）

限界費用とマークアップを次のように生成し、価格を計算せよ。

$$
mc_{jt} = 0.5 + 0.8 w_{jt} + \text{誤差}_{jt}, \quad \text{誤差}_{jt} \sim N(0, 0.3^2)
$$

$$
\text{markup}_{jt} = \max(1.0 + 1.3\, \xi_{jt},\ 0.05)
$$

$$
p_{jt} = mc_{jt} + \text{markup}_{jt}
$$

（`pmax()`を使うと、マークアップが負にならないようにできる。）

生成した価格$p$と需要ショック$\xi$の相関係数を計算し、正の値になっていることを確認せよ。

```{r}
#| label: q2-skeleton
#| eval: false
mc_mkt <- ___ + ___ * w_mkt + matrix(rnorm(___, sd = ___), nrow = ___, ncol = ___)
markup_mkt <- pmax(___, 0.05)
price_mkt <- ___ + ___

cor(as.vector(price_mkt), as.vector(xi_mkt))
```

::: {.callout-tip}
## ヒント：なぜこの生成方法で内生性が生まれるのか

markupが$\xi$に正に反応するように作っているので、$\xi$が高い（消費者に人気がある）財ほど価格も高くなる。これは講義で説明した「企業は$\xi$を知って強気の価格をつける」という状況をそのままコードにしたものである。もしここで`markup_mkt <- 1.0`（定数）としてしまうと、価格は$w$だけで決まり、$\xi$と無相関になる。これが「内生性のない」ベンチマークになる。余裕があれば、後で`markup_mkt`の$\xi$への反応係数を0にした場合にOLSとIVがどう変わるか、比較してみるとよい。
:::

### Q3（10点）：nested logitシェアの計算とデータフレーム化

次のnested logitシェア計算関数を使い（講義ノートの`nested_logit_shares()`と同じもの。コピーして使ってよい）、市場$\times$財のシェアを計算せよ。ネスト構造は`nests <- list(c(1,2), c(3,4))`（財1・2が甘め系ネスト、財3・4が無糖・すっきり系ネスト）とする。

```{r}
#| label: nested-logit-fn-given
nested_logit_shares_batch <- function(x, price, xi, alpha, beta, rho, nests) {
  # x, price, xi: T x J 行列
  delta <- beta * x - alpha * price + xi
  T_ <- nrow(delta); J_ <- ncol(delta)
  lambda <- 1 - rho
  n_nests <- length(nests)

  IV <- matrix(NA, T_, n_nests)
  for (g in seq_len(n_nests)) {
    idx <- nests[[g]]
    IV[, g] <- rowSums(exp(delta[, idx, drop = FALSE] / lambda))
  }
  denom <- 1 + rowSums(IV^lambda)
  s0 <- 1 / denom
  s_g <- (IV^lambda) / denom

  s <- matrix(NA, T_, J_)
  s_within <- matrix(NA, T_, J_)
  for (g in seq_len(n_nests)) {
    idx <- nests[[g]]
    for (j in idx) {
      s_within[, j] <- exp(delta[, j] / lambda) / IV[, g]
      s[, j] <- s_within[, j] * s_g[, g]
    }
  }
  list(s = s, s0 = s0, s_within = s_within)
}
```

この関数を使い、以下を実行せよ。価格内生性とnested-logitの仕様差を混同しないため、**同じ $x,p,\xi$ から2つのDGPを分けて作る**。

1. 真値`alpha=2.0, beta=1.5, rho=0.5`で**nested DGP**のシェアを計算する。
2. 同じ変数を使い、`rho=0`で**plain-logit DGP**のシェアも計算する。
3. nested DGPを`market_df`、plain DGPを`market_df_plain`として、どちらも1行1観測（$T\times J=400$行）にまとめる。
4. `y_berry = log(share) - log(share0)`を両方に、`ln_swithin = log(s_within)`をnested DGPに追加する。

```{r}
#| label: q3-skeleton
#| eval: false
nests <- list(c(1, 2), c(3, 4))
alpha_true <- 2.0
beta_true <- 1.5
rho_true <- 0.5

nl_result <- nested_logit_shares_batch(x_mkt, price_mkt, xi_mkt, alpha_true, beta_true, rho_true, nests)
pl_result <- nested_logit_shares_batch(x_mkt, price_mkt, xi_mkt, alpha_true, beta_true, 0, nests)

market_df <- data.frame(
  market = rep(1:T_mkt, times = J_mkt),
  product = rep(1:J_mkt, each = T_mkt),
  x = as.vector(x_mkt),
  w = as.vector(w_mkt),
  xi = as.vector(xi_mkt),
  price = as.vector(price_mkt),
  share = as.vector(___),
  share0 = rep(___, times = J_mkt),
  s_within = as.vector(___)
) %>%
  mutate(y_berry = ___, ln_swithin = ___)

# plain-logit DGP（同じ x, w, xi, price。違うのは rho=0 のシェア写像だけ）
market_df_plain <- market_df %>%
  mutate(
    share = as.vector(___),
    share0 = rep(___, times = J_mkt),
    y_berry = log(share) - log(share0)
  )
```

### Q4（5点）：記述統計と可視化

1. 財ごとの平均シェア・平均価格を表にまとめよ（`group_by(product) %>% summarize(...)`）。
2. 価格とシェアの散布図を、財ごとに色分けして描け。
3. outside optionの平均シェアを計算せよ。

::: {.callout-tip}
## ヒント

平均シェアはだいたい5〜10%程度、outsideのシェアは60%前後になっているはずである。極端に0%や100%に近い値になっている場合は、Q1・Q2のパラメータ（特に`x_mkt`の平均）を疑ってみるとよい。
:::

## Part 2：推定・推論（60点）

### Q5（10点）：plain logitのBerry反転をOLSで推定する

まず、Q3で別に作った**plain-logit DGP**を使い、Berry反転をOLSで推定する。ここではモデルの関数形は正しいので、真値との差は価格内生性に由来する。nested DGPをplain logitとして推定すると、価格内生性と仕様誤り（within-share項の欠落）が混ざってしまうので、この設問では使わない。

$$
y_{jt} = \ln s_{jt} - \ln s_{0t} = x_{jt}'\beta - \alpha p_{jt} + \xi_{jt}
$$

`fixest::feols()`を使ってOLS推定し、$\hat\alpha$を真値と比較せよ。

```{r}
#| label: q5-skeleton
#| eval: false
ols_plain <- feols(y_berry ~ price + x, data = market_df_plain, vcov = "hetero")
summary(ols_plain)

alpha_ols_plain <- -coef(ols_plain)["price"]
```

### Q6（10点）：nested logit版のBerry反転をOLSで推定する

nested logit版のBerry反転式

$$
y_{jt} = x_{jt}'\beta - \alpha p_{jt} + \rho \ln \bar{s}_{j|g,t} + \xi_{jt}
$$

をOLSで推定せよ（`ln_swithin`を説明変数に追加するだけでよい）。$\hat\alpha$と$\hat\rho$を真値（$\alpha=2.0$, $\rho=0.5$）と比較し、Q5のplain logitの結果と合わせて、**どちらのパラメータがどの方向にどれだけ歪むか**を記述せよ。

```{r}
#| label: q6-skeleton
#| eval: false
ols_nested <- feols(y_berry ~ price + ln_swithin + x, data = market_df, vcov = "hetero")
summary(ols_nested)
```

::: {.callout-tip}
## ヒント：なぜwithin shareもOLSで歪むのか

`ln_swithin`（$\ln \bar{s}_{j|g,t}$）は、同じネスト内の財の$\delta$（したがって$\xi$）にも依存して決まる。つまり価格と同様に、この変数も$\xi_{jt}$と相関する内生変数である。したがって$\alpha$だけでなく$\rho$の推定にもバイアスがかかる。バイアスの方向を予想してから、実際にコードを実行して答え合わせをしてみよう。
:::

### Q7（20点）：2SLSで$\alpha$と$\rho$を回収する

Q6のOLS推定には2つの内生変数（`price`と`ln_swithin`）があるので、含まれる外生変数$x$を除いた操作変数行列に、少なくとも2方向の独立な変動が必要である。この課題では次の2つを除外操作変数として使う。

1. `w`：コストシフター（主に`price`を動かす）
2. `rx_within`：**同じネスト内の他財の商品属性の合計**（相対効用を通じて主に`ln_swithin`を動かす）

`rx_within`は、財$j$がネスト$g$に属するとき、次のように計算する。

$$
\text{rx\_within}_{jt} = \sum_{k \in g,\, k \neq j} x_{kt}.
$$

これは、他製品属性を使うという意味でBLP型IVに着想を得た、この課題専用の操作変数である。財$j$自身の$\xi_{jt}$には影響しない一方、同じネストの他財の魅力度を変え、財$j$との**ネスト内の相対効用**を通じて$\bar s_{j|g,t}$を動かす。

一方、異なるネストの属性和`rx_other`はこの課題では使わない。canonical nested logitでは

$$
\bar s_{j|g,t}
=\frac{\exp(\delta_{jt}/(1-\rho))}
{\sum_{k\in g}\exp(\delta_{kt}/(1-\rho))}
$$

なので、条件付きシェアは他ネストの属性に依存しない。他ネストの属性はネスト選択確率$s_{g,t}$を動かしうるが、$\bar s_{j|g,t}$は動かさないため、このDGPで`ln_swithin`用の除外IVにはならない。標準的なBLP IVは企業IDを使って「同じ企業の他製品」と「他企業の製品」の属性和を分けるが、この課題のネスト区分は企業区分ではないので、両者を同一視しないこと。

```{r}
#| label: q7a-skeleton
#| eval: false
# 財ごとの「同ネスト内の競合のx合計」を計算する関数
compute_within_rival_x <- function(x, nests) {
  T_ <- nrow(x); J_ <- ncol(x)
  nest_of <- integer(J_)
  for (g in seq_along(nests)) {
    for (j in nests[[g]]) nest_of[j] <- g
  }
  rx_within <- matrix(NA, T_, J_)
  for (j in 1:J_) {
    g <- nest_of[j]
    rivals_within <- setdiff(nests[[g]], j)
    rx_within[, j] <- if (length(rivals_within) > 0) rowSums(x[, rivals_within, drop = FALSE]) else 0
  }
  rx_within
}

rx_within_mkt <- compute_within_rival_x(x_mkt, nests)
market_df$rx_within <- as.vector(___)
```

`fixest::feols()`で複数内生変数の2SLSを実行するには、`|`の後ろに`(内生変数1 + 内生変数2) ~ 操作変数1 + 操作変数2 + ...`の形式で書く。

```{r}
#| label: q7b-skeleton
#| eval: false
iv_nested <- feols(
  y_berry ~ x | price + ln_swithin ~ w + rx_within,
  data = market_df, vcov = "hetero"
)
summary(iv_nested)
```

$\hat\alpha$、$\hat\rho$を真値と比較せよ。さらに、含まれる外生変数$x$を条件づけた2本のfirst stageを推定し、(i) 除外IV（`w`, `rx_within`）の係数から作る$2\times2$行列のrankが2であること、(ii) 各first stageにおける2本の除外IVのjoint partial Fを確認せよ。回帰全体のFを報告しないこと。rankが2であることは母集団のrank conditionに対応する標本上の確認であり、partial Fが十分かどうかは弱識別を点検するための別の問いである。

::: {.callout-tip}
## ヒント：係数の取り出し方

`coef(iv_nested)`を実行すると係数名が確認できる。`fixest`のバージョンによって、内生変数の係数名が元の変数名（`price`）のままの場合と、接頭辞が付く場合がありうる。`names(coef(iv_nested))`をまず`print`して名前を確認し、それに応じて`coef(iv_nested)["price"]`または該当する名前でアクセスするとよい。次のような、名前の部分一致で探す頑健なヘルパー関数を用意しておくと安全である。

```{r}
#| label: coef-helper
#| eval: false
get_coef_for <- function(model, varname) {
  cf <- coef(model)
  hit <- grepl(varname, names(cf), fixed = TRUE)
  cf[hit][1]
}
```
:::

::: {.callout-tip}
## ヒント：first-stageの見方

確実な計算方法は、各内生変数について、reduced model（`~ x`）とfull model（`~ x + w + rx_within`）を`lm()`で推定し、`anova(reduced, full)`の2行目のFを読むこと。この古典的partial Fは「2本の除外IVを同時に落とせるか」を検定する。さらに2本のfull first stageから、行を内生変数（`price`, `ln_swithin`）、列を除外IV（`w`, `rx_within`）とする係数行列を作り、`qr(B)$rank`を確認する。`fitstat(fs, "f")`が返す回帰全体のFは、強い$x$だけで大きくなりうるため代用にならない。なお、このrank確認と古典的partial Fは入門的な診断である。異分散と複数内生変数を本格的に扱う実証では、Sanderson–Windmeijer conditional FやKleibergen–Paap統計量、弱IV頑健推論を使う。
:::

::: {.callout-tip}
## ヒント：このDGPで各操作変数は何を動かすか

この課題の価格式にはライバル属性を入れていないので、`rx_within`は構築上、主に`ln_swithin`を動かす。`price`の関連性は主に`w`が担う。ただし`w`が自財価格を通じてwithin shareも動かすため、単に「IVを1本ずつ内生変数へ割り当てる」のではない。2本のfirst stageを同時に見て、除外IV係数行列がfull rankかを確認する必要がある。操作変数の本数が内生変数以上という本数の確認だけでは、識別も強さも保証されない。
:::

### Q8（10点）：弾力性行列の比較（plain/nested × OLS/IV）

次の4通りの弾力性行列を計算し、比較せよ。

1. plain-logit DGP（`market_df_plain`）をOLS推定したモデル
2. plain-logit DGP（`market_df_plain`）を`w`でIV推定したモデル
3. nested-logit DGP（`market_df`）をOLS推定したモデル
4. nested-logit DGP（`market_df`）を`w + rx_within`でIV推定したモデル

plain logitのIVモデルを`market_df_plain`で追加推定せよ。4モデルとも、価格と属性の評価点には各DGPに共通する財別平均を使い、**各推定モデル自身の**$\hat\beta,\hat\alpha,\hat\rho$（plainでは$\rho=0$）から基準確率を計算して数値微分すること。nestedモデルの基準確率に真の$\beta$や観測されたnestedシェアを流用してはいけない。表または図でまとめ、「内生性を無視したときに弾力性と代替性をどう見誤るか」を2〜3行で述べよ。

```{r}
#| label: q8-skeleton
#| eval: false
iv_plain <- feols(y_berry ~ x | price ~ w, data = market_df_plain, vcov = "hetero")

mean_inputs <- market_df %>%
  group_by(product) %>%
  summarize(x = mean(x), price = mean(price), .groups = "drop")

# nested_logit_shares_batch()で各モデル自身の基準確率を作り、価格を1財ずつ微小に動かす
# 4モデルの (beta_hat, alpha_hat, rho_hat) を混同しないこと
```

### Q9（5点）：$\rho = 0$のWald検定

2SLSで得た$\hat\rho$とその標準誤差を使い、$H_0: \rho = 0$（「nested構造は不要でplain logitで十分」という帰無仮説）に対するWald検定統計量を計算せよ。ここでは線形IV式を$\rho<0$にも拡張した無制約推定の通常のWald近似を使う。構造制約$0\leq\rho<1$を推定時に課す場合、$\rho=0$は境界なので通常の$\chi^2_1$近似には注意が要る。

$$
W = \left(\frac{\hat\rho}{\text{se}(\hat\rho)}\right)^2
$$

は自由度1のカイ二乗分布に従う。5%水準の臨界値（3.84）と比較し、$H_0$を棄却できるか判定せよ。

```{r}
#| label: q9-skeleton
#| eval: false
get_se_for <- function(model, varname) {
  s <- se(model)
  hit <- grepl(varname, names(s), fixed = TRUE)
  s[hit][1]
}

rho_hat <- get_coef_for(iv_nested, "ln_swithin")
se_rho <- get_se_for(iv_nested, "ln_swithin")
wald_stat <- (rho_hat / se_rho)^2
```

### Q10（5点）：経営への含意（3行で）

あなたは飲料メーカーの商品企画担当として、この推定結果をもとに次の問いに答えなければならない。「甘め系ネストの財1を値上げしたとき、一番シェアを奪われる（＝一番怖い競合になる）のはどの財か」。2SLSで推定した弾力性行列を根拠に、3行以内で答えよ。OLSに基づく弾力性行列と比べ、競合順位が変わるかだけでなく、自己反応・同ネスト内への流出・別ネストへの流出の大きさがどう違うかにも触れること。

::: {.callout-warning}
## 提出前のチェック

- チャンクを上から順に実行して、未定義変数によるエラーが出ないか確認したか。
- `alpha_ols_plain`、`alpha_ols_nested`（Q6の$\hat\alpha$）、`alpha_iv`（Q7の$\hat\alpha$）を明確に区別して報告しているか。
- Q8の弾力性行列は、対角成分（自己弾力性）と非対角成分（交差弾力性）の両方に触れているか。
:::
