課題2：動画配信サービスのプラン選択とlogit需要推定

計量経済学II

作者

Kei Ikegami

重要この課題について

この課題は、講義ノートで扱った「缶コーヒー棚」のロジックを、動画配信サービスのプラン選択という別の題材で自分の手で再現する。属性、パラメータの値、価格帯を変えてあるので、DGPと推定の中身を理解しながら実装してほしい。講義の特性アプローチと同じく、全員・全プランに共通する属性の係数から各プランの効用を作り、選択確率を需要関数として読む。

配点は合計100点。Part 1はDGPの実装と可視化（30点）、Part 2は推定・推論・反実仮想（70点）。手を動かして2〜4時間程度で終わる分量を想定している。ヒントは惜しみなく置いてあるので、詰まったら読んでほしい。

設定：動画配信サービスのプラン選択

とある動画配信サービスには、3つの有料プラン（Basic・Standard・Premium）と、「加入しない」（outside option）という選択肢がある。各プランには、月額料金（price、円）、画質（quality、1=SD、2=HD、3=4Kの3段階）、同時視聴可能数（screens、同時に視聴できる画面の数）という3つの属性がある。

消費者 \(i\) は、これらのプランのうちどれに加入するか（あるいは加入しないか）を、講義で学んだのと同じrandom utility modelで決める。消費者 \(i\) がプラン \(j\) から得る効用は

\[ u_{ij} = v_{ij} + \varepsilon_{ij}, \qquad v_{ij} = \beta_{\text{quality}} \cdot \text{quality}_j + \beta_{\text{screens}} \cdot \text{screens}_j - \alpha \cdot \text{price}_{ij} \]

である。\(\beta_{\text{quality}}\)、\(\beta_{\text{screens}}\)、\(\alpha\) は全員・全プランに共通とし、\(\varepsilon_{ij}\) は消費者間・選択肢間で独立なスケール1のGumbel分布（Type I極値分布）に従う。提示価格の変動も、この誤差と独立に生成する。加入しない（outside, \(j=0\)）の決定的効用を \(v_{i0}=0\) とする位置正規化に加え、この誤差スケールの固定が尺度正規化である。outsideの総効用は \(u_{i0}=\varepsilon_{i0}\) であり、0に固定するわけではない。推定する係数は、誤差のスケールに対する相対値として解釈する。

同じ価格ベクトル \(p=(p_{\mathrm{Basic}},p_{\mathrm{Standard}},p_{\mathrm{Premium}})\) に直面する市場では、logitシェアを \(s_j(p)\)、市場規模を \(M\) とすると、各プランの需要は \(q_j(p)=Ms_j(p)\) である。自分の価格だけでなく、他のプランの価格も分母を通じて需要に影響する。

なお、このDGPにはブランド固有の切片（講義ノートのconst_A, const_Bのようなプラン別の定数項）を置いていない。理由は設問1.4で考えてもらう。

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

1.1 プラン属性とパラメータの設定（5点）

次のプラン属性と真のパラメータ値を使う。

N <- 3000

# プラン属性: Basic, Standard, Premium
price_base <- c(700, 1200, 1800)   # 月額（円）
quality    <- c(1, 2, 3)            # 1=SD, 2=HD, 3=4K
screens    <- c(1, 2, 4)            # 同時視聴可能数

true_params <- list(
  beta_quality = 0.9,
  beta_screens = 0.35,
  alpha_price  = 0.003   # 価格は円単位。0.003 x 1200円 ~ 3.6程度の効用差になるスケール
)

設問1.1：alpha_price = 0.003という値が、月額700〜1800円という価格帯に対してどの程度のスケール感か、alpha_price * price_baseを計算して確認せよ。もしalpha_price = 0.03だったら（10倍大きい場合）、効用差はどうなるか、シェアにどう影響しそうか、1〜2文で予想してから次に進むこと。

# ヒント: 以下を実行してスケール感を確認する
true_params$alpha_price * price_base

1.2 Gumbel乱数によるデータ生成（15点）

講義で導入した通り、Gumbel分布に従う乱数は一様乱数から-log(-log(runif(n)))で作れる。各消費者について、わずかな価格変動（プロモーションなど）を加えたうえで、効用を計算し、最大効用のプランを選ばせる関数を書け。

draw_gumbel <- function(n) {
  u <- runif(n)
  -log(-log(u))
}

simulate_subscription_choice <- function(n, price_base, quality, screens, params) {
  # ヒント1: 価格にわずかな個人間ばらつきを与える（プロモーションや地域差のイメージ）。
  #   例: price_noise <- matrix(rnorm(n * 3, mean = 0, sd = 40), nrow = n, ncol = 3)
  #       price <- sweep(price_noise, 2, price_base, "+")
  #       price <- pmin(pmax(price, 400), 2500)  # 極端な値をクリップ

  # ヒント2: 決定的効用 v_ij を3プラン分計算する（cbindで3列の行列にする）。
  #   delta <- cbind(
  #     params$beta_quality*quality[1] + params$beta_screens*screens[1] - params$alpha_price*price[,1],
  #     ...
  #   )

  # ヒント3: outside option（決定的効用0）の列を先頭に追加する。
  #   delta_full <- cbind(0, delta)

  # ヒント4: Gumbel誤差を加えて、行ごとに最大値のインデックスを選ばせる。
  #   eps <- matrix(draw_gumbel(n * 4), nrow = n, ncol = 4)
  #   u <- delta_full + eps
  #   choice <- apply(u, 1, which.max) - 1   # 0=outside, 1=Basic, 2=Standard, 3=Premium

  # 戻り値として choice と price の両方を返すこと（後で使う）
}
ヒントヒント：DGPのデバッグ手順
  1. まずn=10くらいの小さいサンプルでコードを走らせ、choiceベクトルの中身が0〜3の整数になっているか確認する。
  2. table(choice)で頻度を見て、極端に偏っていないか（例えば全員が0=outsideになっていないか）を確認する。全員が同じ選択肢に集中する場合、たいてい価格や属性のスケールがおかしい（1.1の設問を振り返ること）。
  3. うまくいったらn=3000で本番実行する。

シミュレーションしたデータで、選択シェアの棒グラフを描け。

sim_data <- simulate_subscription_choice(N, price_base, quality, screens, true_params)
choice_vec <- sim_data$choice
price_mat <- sim_data$price

# ヒント: table()とggplotのgeom_col()で、outside/Basic/Standard/Premiumのシェアを棒グラフにする

1.3 属性別クロス集計（5点）

各プランの選択者数（人数とシェア）を表にまとめよ。さらに、価格帯別（例えば価格を安い/普通/高いの3群に分けるなど、自分で工夫してよい）に選択シェアがどう変わるかを簡単なクロス集計で確認せよ。

ヒントヒント：クロス集計の作り方

dplyr::group_by()とsummarise()を使うのが簡単である。価格変動が小さいシミュレーション設計であれば、価格帯の違いはそれほど劇的にシェアを変えないかもしれない。それならそれで「なぜそうなるか」を一言添えるとよい（プロモーション幅が小さいから、など）。

1.4 考察：なぜプラン固有の切片を入れなかったか（5点）

このDGPでは、3プランの非価格部分を2つの属性係数だけで説明する。講義の最初の推定例にあったconst_A, const_Bのように、BasicとStandardの切片を両方追加すると、なぜ属性係数と分けて推定できなくなるか。2〜3文で説明せよ。ヒント：qualityとscreensは各プランに固定されている。outside/Basic/Standard/Premiumを行、Basicダミー・Standardダミー・quality・screensを列にした行列のランクと列数を比べる。これは「切片を1つでも入れると必ず識別できない」という主張ではない。切片を置かない現モデルは、属性で説明できない恒常的なプランの魅力はない、という制約でもある。

Part 2：推定・推論（70点）

2.1 自作対数尤度とoptim()（20点）

講義ノートのneg_loglik_mnl関数を参考に、このプラン選択データに対する対数尤度関数を自作し、optim(method="BFGS", hessian=TRUE)でパラメータ\((\beta_{\text{quality}}, \beta_{\text{screens}}, \alpha)\)を推定せよ。真値と推定値を比較する表を作ること。動作確認として、同じ係数の候補を使い、最初の3人が実際に選んだ選択肢の対数確率を個別に足した値と、関数が返す負の対数尤度の符号を反転した値が一致することを確認せよ。

build_design_sub <- function(price, quality, screens) {
  # ヒント: 講義ノートのbuild_design()と同じ発想。
  # ただしこの課題では選択肢固有の切片が無いので、パラメータは3つだけ
  # （beta_quality, beta_screens, alpha）。
  # X[, j, ] の3列がそれぞれ (quality_j, screens_j, -price_j) になるように配列を組む。
  n <- nrow(price)
  X <- array(0, dim = c(n, 4, 3))
  # 第2次元の1列目 = outside（全て0）、2列目 = Basic、3列目 = Standard、4列目 = Premium
  # X[, 2, 1] <- quality[1]; X[, 2, 2] <- screens[1]; X[, 2, 3] <- -price[, 1]
  # ... (Standard, Premiumも同様に埋める)
  X
}

neg_loglik_sub <- function(theta, X, choice) {
  # ヒント: 講義ノートのneg_loglik_mnl()とほぼ同じ構造。
  # V <- 各選択肢のX[,j,] %*% theta を並べた n x 4 行列
  # オーバーフロー対策で各行の最大値を引く
  # 選択された選択肢の対数確率の合計にマイナスをつけて返す
}
ヒントヒント：真値回収がうまくいかないとき

推定値が真値から大きくズレる場合、まず疑うべきは次の3点である。

  1. build_design_subの列の並びと、thetaの要素の並びが対応しているか（順番がズレていないか）。
  2. choiceが0始まり（0=outside）になっているか、それとも1始まりになっているか。Rのインデックスは1始まりなので、choice + 1のような変換を忘れていないか。
  3. optimの初期値が極端すぎないか（c(0,0,0)のようなニュートラルな値から始めるのが無難）。

価格係数は他の係数より桁が小さい。この課題ではcontrol = list(parscale = c(1, 1, 0.001), ndeps = c(0.001, 0.001, 0.00001), maxit = 1000, reltol = 1e-10)を指定する。parscaleで探索の尺度をそろえ、ndepsで数値微分の幅を調整する。価格が円単位で大きいため、微分幅が粗いと係数が近くてもHessianからの標準誤差に誤差が出る。opt_result$convergence == 0も確認すること。

2.2 HessianからのSE・信頼区間・検定（15点）

optimが返すhessianを使って、標準誤差と95%信頼区間を計算せよ。また、「価格は需要に影響しない」という帰無仮説 \(H_0: \alpha = 0\) をWald検定せよ（大標本の\(z\)統計量を計算し、有意かどうか判断する）。ここでの逆Hessianは、正しく特定されたiid選択尤度のもとでのモデルベースSEであることも1文で明記せよ。

ヒントヒント：仮説検定の作り方

大標本のWald \(z\)統計量は \(\hat\alpha / \text{SE}(\hat\alpha)\) で計算できる。\(|z|>1.96\) なら、標準正規近似による両側5%水準で \(H_0\) を棄却する。

2.3 mlogitパッケージでの検算（10点）

mlogitパッケージを使って同じモデルを推定し、自作の推定値と近い値になることを確認せよ。提出用レンダーでは#| eval: falseで構わないが、可能なら手元のRセッションで一度実行すること。eval: falseのチャンクはレンダー時には実行されない。

library(mlogit)

# ヒント: long形式に変換し、mlogit.data(..., chid.var = "id")で各人の選択集合を指定する。
# 各idについてchosenが4行中ちょうど1つだけTRUEか、推定前に必ず確認する。
# 講義ノートのlec2-mlogit-checkチャンクの書き方を参考にすること。

2.4 反実仮想：Standardプランを値上げしたら（15点）

この設問では、全員に共通の基準価格 price_base = c(700, 1200, 1800) を提示する市場を考える。Standardプランの月額だけを+200円した場合の、各プランの予測シェアを計算せよ。市場規模を \(M=10{,}000\) 人とし、値上げ前後のシェアと予測加入者数を表にする。さらに、Standardの予測需要の減少分に対して、Basic・Premium・outsideのそれぞれの増加分が占める割合を計算せよ。outsideは「加入しない」であり、この1回の選択データだけから実際の解約者を特定できるわけではない。

値上げ前の価格・シェアを使い、Standardの自己価格弾力性 \(\eta_{SS}=-\hat\alpha p_S(1-s_S)\) と、Standardの価格に対するBasicの交差価格弾力性 \(\eta_{BS}=\hat\alpha p_Ss_S\) を求め、Standardを1%値上げしたときの需要の変化をそれぞれ1文で解釈せよ。ここでの微分による弾力性は小さな価格変化に対する近似であり、+200円の予測はlogit式に新価格を代入して計算する。

その上で、この「流出の内訳」に何か規則性がないかを観察し、1〜2文で言葉にしてみよ（ヒント：それぞれの流出先の、値上げ前のシェアの大きさと比べてみること）。これは、講義で紹介したIIA（Independence of Irrelevant Alternatives）という性質の確認である。次回は、この性質が代替パターンに課す制約を詳しく考える。

ヒントヒント：反実仮想の計算手順
  1. price_ref <- matrix(price_base, nrow = 1)を作り、この共通価格でのシェアを計算する（推定したtheta_hatを使う）。ここでは個人別価格の予測を平均しない。
  2. Standardの価格だけ+200円した新しい価格行列を作る。
  3. 新しい価格でのシェアを計算する。
  4. 2つのシェアの差を取る。Standardのシェアがどれだけ減ったか、他のプランのシェアがどれだけ増えたかを見る。
  5. 「Standardからの流出分」に対する「各流出先が受け取った割合」を redistribution_share <- gain_others / lost_from_standard で計算する。あわせて、値上げ前の各行き先のシェアをStandard以外のシェア合計で割った条件付きシェアを作り、2つを同じ表で比べる。

共通の価格・属性に条件付けたlogitでは、この2つの割合が一致する。個人ごとに異なる価格で予測し、そのシェアを平均した場合には、集計後も厳密に一致するとは限らない。

2.5 モンテカルロによる被覆率の確認（5点）

同じDGPから200回データを生成し直し、そのたびに推定・95%信頼区間の計算を行い、真値が信頼区間に入る頻度（被覆率）を計算せよ。被覆率が0.95に近いかどうかをコメントすること。計算が収束しなかった反復や、有効な標準誤差を得られなかった反復があれば、その件数と被覆率の計算に使った件数を報告せよ。

ヒントヒント：モンテカルロループの書き方

講義ノートのlec2-montecarloチャンクとほぼ同じ構造でよい。ループの中でsimulate_subscription_choice()とoptim()を繰り返し呼び出し、都度信頼区間を計算して、真値が入っているかをTRUE/FALSEで記録し、最後に平均を取ればよい。200回のループは数十秒程度かかることがあるので、気長に待つこと。

2.6 経営への含意（5点）

以上の分析結果を踏まえて、動画配信サービスの経営陣に向けて、「Standardプランの値上げを検討する際に注意すべきこと」を3行程度で書け。価格弾力性、カニバリゼーション（自社プラン間での顧客の奪い合い）、outside（非加入）の増加を踏まえること。講義の1商品だけを売る企業と違い、ここでは同じ企業が3プランを販売している。どの価格が利益を最大にするかは、この需要予測に加えて各プランの費用も使って判断する必要がある。

  • Part 1（30点）：1.1 スケールの確認（5点）、1.2 DGPとシェアの棒グラフ（15点）、1.3 クロス集計（5点）、1.4 識別の考察（5点）。
  • Part 2（70点）：2.1 自作MLE・3人分の尤度確認・真値回収表（20点）、2.2 SE・CI・Wald検定（15点）、2.3 mlogitでの検算（10点）、2.4 反実仮想・弾力性・IIAの観察（15点）、2.5 モンテカルロ被覆率（5点）、2.6 経営含意3行（5点）。
  • 部分点の基準：DGPのコードが動かなくても、ロジックが正しく記述されていれば部分点を与える。真値回収がやや外れていても（サンプリング誤差の範囲内であれば）、実装の考え方が正しければ大きく減点しない。

提出方法

このqmdファイルに直接コードと回答を書き込み、レンダリングしたHTMLファイルと合わせて提出すること。