課題3：nested logitの市場レベル需要推定と価格内生性

計量経済学II

作者

Kei Ikegami

重要この課題のゴール

第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チャンクで既に固定されているので、以降このチャンクの実行順序を変えないこと。

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 = ___)
ヒントヒント

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\)の相関係数を計算し、正の値になっていることを確認せよ。

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))
ヒントヒント：なぜこの生成方法で内生性が生まれるのか

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が無糖・すっきり系ネスト）とする。

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に追加する。
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の平均シェアを計算せよ。
ヒントヒント

平均シェアはだいたい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\)を真値と比較せよ。

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の結果と合わせて、どちらのパラメータがどの方向にどれだけ歪むかを記述せよ。

ols_nested <- feols(y_berry ~ price + ln_swithin + x, data = market_df, vcov = "hetero")
summary(ols_nested)
ヒントヒント：なぜ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を使って「同じ企業の他製品」と「他企業の製品」の属性和を分けるが、この課題のネスト区分は企業区分ではないので、両者を同一視しないこと。

# 財ごとの「同ネスト内の競合の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 + ...の形式で書く。

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が十分かどうかは弱識別を点検するための別の問いである。

ヒントヒント：係数の取り出し方

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

get_coef_for <- function(model, varname) {
  cf <- coef(model)
  hit <- grepl(varname, names(cf), fixed = TRUE)
  cf[hit][1]
}
ヒントヒント：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頑健推論を使う。

ヒントヒント：この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行で述べよ。

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\)を棄却できるか判定せよ。

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に基づく弾力性行列と比べ、競合順位が変わるかだけでなく、自己反応・同ネスト内への流出・別ネストへの流出の大きさがどう違うかにも触れること。

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