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

計量経済学II

作者

Kei Ikegami

重要この課題について

この課題は、講義（lecture2.qmd）で扱った「缶コーヒー棚」のロジックを、動画配信サービスのプラン選択という別の題材で自分の手で再現する。設定（選択肢の数、パラメータの値、文脈）は講義ノートとは変えてあるので、コピー＆ペーストではなく、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} \]

であり、\(\varepsilon_{ij}\) はスケール1のGumbel分布（Type I極値分布）に従う。加入しない（outside, \(j=0\)）の効用を \(v_{i0}=0\) とする位置正規化に加え、この誤差スケールの固定が尺度正規化である。したがって推定する係数は、誤差のスケールに対する相対値として解釈する。

なお、この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 属性別クロス集計（10点）

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

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

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

1.4 考察：なぜプラン固有の切片を入れなかったか（考察のみ、点数はPart 2の設問と合算）

このDGPには、講義ノートの缶コーヒーの例にあったconst_A, const_Bのような、プラン固有の切片（ダミー変数）を入れていない。なぜだと思うか、2〜3文で考えてみよ。ヒント：この課題ではqualityとscreensが各プランに固定されている。outside/Basic/Standard/Premiumを行、候補の切片ダミー・quality・screensを列にした小さな行列を書き、列ランクと列数を比べるとよい。なお、属性が市場や選択機会をまたいで変化するデータでは結論が変わり得る。

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

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

講義ノートのneg_loglik_mnl関数を参考に、このプラン選択データに対する対数尤度関数を自作し、optim(method="BFGS", hessian=TRUE)でパラメータ\((\beta_{\text{quality}}, \beta_{\text{screens}}, \alpha)\)を推定せよ。真値と推定値を比較する表を作ること。

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))
  # 選択肢1 = Basic, 2 = Standard, 3 = Premium（列0=outsideは全て0のまま）
  # 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)のようなニュートラルな値から始めるのが無難）。

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点）

Standardプランの月額を+200円した場合の、各プランの予測シェアを計算せよ。値上げ前後のシェアを比較する表を作り、「値上げによってStandardから離脱した客のうち、どれくらいの割合がBasicに流れ、どれくらいがPremiumに流れ、どれくらいがoutside（解約）に流れたか」を計算せよ。

その上で、この「流出の内訳」に何か規則性がないかを観察し、1〜2文で言葉にしてみよ（ヒント：それぞれの流出先の、値上げ前のシェアの大きさと比べてみること）。この観察は、次回の講義で学ぶIIA（Independence of Irrelevant Alternatives）という性質の伏線になっている。

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

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

同じDGPから200回データを生成し直し、そのたびに推定・95%信頼区間の計算を行い、真値が信頼区間に入る頻度（被覆率）を計算せよ。被覆率が0.95に近いかどうかをコメントすること。

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

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

2.6 経営への含意（合計点に含む）

以上の分析結果を踏まえて、動画配信サービスの経営陣に向けて、「Standardプランの値上げを検討する際に注意すべきこと」を3行程度で書け。価格弾力性の大きさ、カニバリゼーション（自社プラン間での顧客の奪い合い）、outsideへの流出（純粋な解約）の3つの観点を意識すること。

  • Part 1（30点）：DGPが正しく実装されている（15点）、可視化・記述統計が適切（10点）、1.4の考察に触れている（5点）。
  • Part 2（70点）：自作MLEの実装・真値回収表（20点）、SE・CI・Wald検定（15点）、mlogitでの検算（10点）、反実仮想とIIAの兆候の指摘（15点）、モンテカルロ被覆率（10点）。2.6の経営含意3行は、Part 2全体の完成度に対する加点要素として評価する。
  • 部分点の基準：DGPのコードが動かなくても、ロジックが正しく記述されていれば部分点を与える。真値回収がやや外れていても（サンプリング誤差の範囲内であれば）、実装の考え方が正しければ大きく減点しない。

提出方法

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