Lecture 2：離散選択モデル1 — logitで需要を考える

計量経済学II

重要この講義で押さえたいこと

• 製品を特徴の束として表し、消費者の選択から「価格を変えると購入数がどう変わるか」という需要関数を組み立てられる。
• Gumbel分布（Type I極値分布）を効用のランダム項に仮定すると、選択確率が驚くほどシンプルな閉形式（logit公式）になる理由を理解する。
• 価格弾力性から「値上げしたら客がどこへ流れるか」を読み取り、企業の利益最大化や2社の価格競争の均衡につなげられる。
• 最尤法（MLE）の考え方を、3人の購入データから出発してlogitモデルまで自分の手で追える。
• optim()で商品の属性に対する共通の係数も推定し、属性変更後の需要を予測する。シミュレーションで真値回収を確認し、Hessianから標準誤差・信頼区間を計算する。

前回とのつながり・今回のゴール

第1回では、実証IOとquantitative marketingが企業や政策の意思決定にどう結びつくかを概観し、予測・因果効果・反実仮想という3つの問いの違いを確認した。そのなかで、「価格を上げたら需要はどう動くか」という反実仮想の問いに答えるには、単純な回帰だけでは力不足だという話をした。

今回からしばらくは、この「反実仮想に答えるための道具」として離散選択モデル（discrete choice model）を学ぶ。第2回である今回は、その最も基本的な形である multinomial logit（MNL、多項ロジット） を扱う。ゴールは3つ。

1. 製品の品質や好みの違いを効用に書き、消費者の選択を集めて需要関数を作れるようになる（RUM）。
2. なぜlogitという特定の形になるのかをGumbel分布の性質から理解し、需要曲線・需要曲面から自己価格弾力性と交差価格弾力性の意味と式をつかむ。
3. 実際にRでシミュレーションデータを作り、自作の最尤法でパラメータを推定する。商品の属性に対する共通の係数を推定して属性変更後の需要を予測し、価格設定にも使う。

次回（第3回）は、今回作るMNLモデルの限界—IIA（独立性仮定）という代替性の制約と、価格が需要ショックと相関する内生性の問題—に進む。

なぜ回帰では足りないか

あなたはコンビニの飲料担当で、缶コーヒー棚を見ている。棚には3つのブランドの缶コーヒーが並んでいる（ブランドA、B、Cと呼ぶことにする）。お客さんは、この3本のうちどれかを買うか、あるいは「今日はやめとくか」と何も買わずに立ち去るか（買わないという選択肢、outside option）のいずれかを選ぶ。

ここでの「アウトカム」は、売上金額でも購入量でもない。「どれを選んだか」という、カテゴリカルな1つの選択である。あなたが知りたいのは、次のような問いだ。

ブランドAを10円値上げしたら、Aの購入数は何%減るか。そして、その減った分の客は、Bに流れるのか、Cに流れるのか、それとも「買わない」を選ぶのか。

「Aの売上」を「Aの価格」に回帰する単品の1本の式だけでは、この問いには答えられない。価格弾力性らしき数字は出ても、その客がどこへ流れたかは式に入っておらず、各選択肢の確率が非負で合計1になることも保証されない。

こんな時に離散選択モデルが便利である。効用最大化から、A・B・C・outsideのシェアが同時にどう動くかを導くことができるからだ。

注意点として、個人のアウトカムが「どれを選ぶか」でも、それを人数分集めると、各商品の需要量になる。例えば、1,000人がそれぞれ最大1本を買う市場で、Aを選ぶ確率が0.4なら、Aの予測購入数は \(1{,}000\times0.4=400\) 本である。Aを値上げしたとき、その確率が0.3になるモデルなら、予測購入数は300本になる。離散選択モデルで推定したいのは、まさにこの「価格から購入確率を経て購入数へ」という関係である。

ノートビジネスの現場で：コンビニの缶コーヒー棚

棚割り（シェルフスペースの配分）や価格改定の意思決定は、まさにこの「シェアの再配分」を予測する仕事である。POSデータには「いつ、どの店で、どの商品が、いくらで売れたか」の記録があるが、そこから「Aを値上げしたときにBへの流出がどれくらいか」を読み取るには、単品ごとの売上回帰ではなく、選択肢全体を条件付き分布として捉えるモデルが必要になる。この講義で作るのは、まさにその最小の道具である。

Random Utility Model（RUM）

計量経済学Iで、「モデルとは条件付き分布の束である」という言い方をした。回帰分析であれば、\(Y \mid X\) の条件付き分布（の平均）を書くのがモデルだった。離散選択モデルでも同じ発想を使う。今回は、アウトカムが「選んだ財」であり、条件は「各財の属性（価格、甘さなど）」になる。

消費者 \(i\) が、財 \(j = 0, 1, \ldots, J\) の中から1つを選ぶとする（\(j=0\) は「何も買わない」というoutside option）。消費者 \(i\) が財 \(j\) を選んだときに得る効用を、次のように分解する。

\[ u_{ij} = v_{ij} + \varepsilon_{ij} \]

• \(v_{ij}\)：分析者が観測できる属性（価格、甘さ、ブランドなど）から決まる、決定的な（deterministic）効用の部分。
• \(\varepsilon_{ij}\)：分析者には観測できない要因による、ランダムな効用の部分。

需要を考えるため、決定的効用には価格を明示的に入れる。例えば、非価格属性のベクトルを \(x_j\)、その評価を \(\beta\) として、

\[ v_{ij}=x_j'\beta-\alpha p_j,\qquad \alpha>0 \]

と書く。\(x_j'\beta\) は「甘さ、品質、ブランドなどから得る魅力」、\(-\alpha p_j\) は「高いほど買いたくなくなる」という部分である。他の商品の価格と属性を固定してAだけを10円値上げすると、Aの効用だけが \(10\alpha\) 下がる。ここでは属性への評価 \(\beta\) と価格への反応 \(\alpha\) を全員共通とし、同じ属性と価格のもとでの人による選択の違いを \(\varepsilon_{ij}\) で表現する。つまり、\(v_{ij}\) と書いてるがこのモデルでは \(v_{j}\) である。

消費者は、自分にとって効用が最大の財を選ぶと仮定する。つまり、

\[ i \text{ が財 } j \text{ を選ぶ} \iff u_{ij} \geq u_{ik} \ \text{for all } k = 0,1,\ldots,J \]

これがrandom utility model（RUM）である。「モデル＝条件付き分布の束」という言葉に戻ると、離散選択モデルとは、商品属性を条件として、消費者がどれを選ぶかという選択の条件付き分布を書くことに他ならない。\(\varepsilon_{ij}\) にどんな分布を仮定するかによって、この条件付き分布の形（＝選択確率の式）が変わる。今回、というかlogit 離散選択モデルはその分布として、Gumbel分布 と呼ばれる特殊な分布を使った時に出てくるモデルである。

\(\varepsilon_{ij}\) は「分析者に見えない要因」と言ったが、具体的に何を指すのだろうか。Manski (1977) に沿った古典的な整理では、ランダム性の源泉として次の4つが挙げられる。

1. 観測されない選択肢属性（unobserved alternative attributes）：パッケージ、香り、混雑など、効用には入るが分析者のデータにない属性。
2. 観測されない選好の違い（unobserved taste variation）：同じ属性でも評価が人や機会によって違うこと。その日の気分や体調もここに含められる。
3. 測定誤差・不完全な情報（measurement error / imperfect information）：属性や認知を正確に測れないこと。
4. 代理変数による近似（proxy or instrumental variables）：本当に効用に入る概念を直接測れず、不完全な代理変数で表すこと。ここでいう古い文献の “instrumental” は、因果推論の除外操作変数というより代理指標の意味である。

したがって \(\varepsilon_{ij}\) は、必ずしも「同じ人が理由なくランダムに選ぶ」という行動仮定ではない。分析者に見えない異質性・属性・測定のずれをまとめた項である。どの源泉を想定するかによって、独立性や分散の仮定が妥当かも変わる。

Lecture 1との接続：製品の特徴から需要を作る

第1回では、需要関数を使って、価格や競争の条件を変えたときに何が起こるかを考えた。今回は、その需要関数の中身を消費者の選択から作る。製品の属性と価格から効用を決め、各人の選択を集計して需要量を得る、という順番である。

製品を特徴の束として表す：特性アプローチ

コーヒーを「商品A」「商品B」という名前だけで区別せず、苦味・甘さなどの特徴の束 \(x_j\) として表そう。消費者が製品の持つ特性から効用を得るという発想は、Lancaster (1966) に代表される特性アプローチ（characteristics approach、特性空間アプローチ）と呼ばれる。

この回では、属性を評価する係数 \(\beta\) と価格係数 \(\alpha\) を全員共通にする。例えば、苦味1単位の係数を1、甘さ1単位の係数を0.5、価格係数を \(\alpha=0.01\)（価格は円単位）とすると、共通部分の効用は

\[ v_j=1\times\text{苦味}_j+0.5\times\text{甘さ}_j-0.01p_j \]

となる。この小例では全員が同じ属性と価格に直面するので、\(v_{ij}\) を \(v_j\) と書いた。

商品	苦味	甘さ	価格	共通部分の効用 \(v_j\)
A	2	0	100円	\(2-0.01\times100=1\)
B	0	2	100円	\(0.5\times2-0.01\times100=0\)

ここで計算した1や0は、効用の単位で測った値であり、円で測った評価額ではない。買わない選択肢については共通部分を \(v_0=0\) とする。

共通の係数でも、人によって選ぶ商品は違う

実際の効用は \(u_{ij}=v_j+\varepsilon_{ij}\) である。共通部分ではAがBより高くても、分析者に見えない好みや事情を含む \(\varepsilon_{ij}\) によって、人ごとの選択は分かれる。仕組みを見るために、普段は見えないこの項を3人分だけ仮に見せよう。各人はA・B・買わない、のうち最も効用の高いものを選ぶ。

消費者	\((\varepsilon_{iA},\varepsilon_{iB},\varepsilon_{i0})\)	A・Bとも100円のとき	Aだけ150円にしたとき
1人目	\((0.4,\ 0.2,\ -0.2)\)	A	A
2人目	\((-0.3,\ 0.4,\ 0)\)	A	B
3人目	\((-0.5,\ 0.8,\ 0.1)\)	B	B

例えば2人目の効用は、値上げ前にはAが \(1-0.3=0.7\)、Bが \(0+0.4=0.4\)、買わない場合が \(0+0=0\) なのでAを選ぶ。属性と \(\varepsilon_{ij}\) を固定したままAだけ150円にすると、Aの共通部分が \(2-0.01\times150=0.5\) に下がり、この人のAの効用は \(0.5-0.3=0.2\) となる。今度はBの0.4の方が大きいので、Bに移る。

3人の選択を集計すると、Aの需要は2本から1本へ、Bの需要は1本から2本へ変わる。属性への係数も価格係数も全員共通のままで、個人の選択と市場の需要を結びつけられた。

実際のデータでは、表の \(\varepsilon_{ij}\) は観測できず、一人ひとりの値をこのように知ることはできない。そこで、その分布を仮定して「Aを選ぶ確率」「Bを選ぶ確率」を計算し、それを集計して需要を予測する。このあと導くlogit公式が、その計算を可能にする。

特製アプローチの偉さ： 製品ごとの需要式を、少数の属性の係数で結びつける

例えば、各製品の需要を他の全製品の価格の一次式で直接書くなら、

\[ q_A=b_A+\gamma_{AA}p_A+\gamma_{AB}p_B+\gamma_{AC}p_C \]

という式がBとCについても必要になる。価格係数をすべて自由にすると、3製品で9個、100製品で10,000個になる。こんなことやってられない。（まぁこういうことをやる需要推定の本格的な研究(Almost Ideal Demand Systemとか)もあるが、それは製品数とか製品カテゴリーが少数の時の話である。）

特性アプローチでは、少数の属性の係数とランダム項の分布への仮定によって、多数の製品の需要を結びつけられる。製品が何千個あろうと、その属性の数は元から決まってるのでなんとかなる、みたいな感じである。

これは多くの需要関数を、消費者の選択という共通の仕組みから作るという発想である。その分、どの商品からどの商品へ客が移るかにも制約が入る。まず第2回では共通の係数によるlogitを学び、その代替パターンの制約を第3回で確かめ、属性への評価そのものを人ごとに変える拡張は第4回で扱う。

Gumbel分布とlogit公式

なぜGumbelなのか

\(\varepsilon_{ij}\) の分布として、Type I 極値分布（Gumbel分布）を仮定するのが、離散選択モデルの標準的な出発点である。Gumbel分布の累積分布関数は

\[ F(\varepsilon) = \exp(-\exp(-\varepsilon)) \]

である。なぜこの分布を選ぶのか。理由は主に3つある。

1. 2つのGumbel確率変数の差が、ロジスティック分布に従う。これによって、2財の選択確率が閉形式（ロジスティック関数）で書ける。
2. 独立なGumbel確率変数の最大値は、再びGumbel分布に従う（max-stability）。これにより、\(J+1\)個の選択肢があっても、選択確率の式がJの数に関わらず綺麗な形にまとまる。
3. 結果として、選択確率が積分を含まない閉形式（closed form）で書ける。他の分布（例えば正規分布、これはprobitモデルになる）を仮定すると、選択確率は多次元積分になり、数値計算が格段に重くなる。

この3つの性質のおかげで、次の公式（multinomial logitの選択確率）が得られる。消費者 \(i\) が財 \(j\) を選ぶ確率は、

\[ s_{ij} = \frac{\exp(v_{ij})}{\sum_{k=0}^{J} \exp(v_{ik})} \]

これがlogit公式である。導出の細部は、以下にまとめておく。今は「そういうものだ」と思って先に進んでも構わない。70%ぐらい理解できたら十分、あとで戻ってくればいい。これを70%ルールと呼ぶ。新しいことを学ぶ際の鉄則である。

ノート導出：logitの選択確率

2財版（\(J=1\)、財0と財1のみ）

消費者は \(u_{i1} = v_{i1} + \varepsilon_{i1}\) と \(u_{i0} = v_{i0} + \varepsilon_{i0}\) を比べて、大きい方を選ぶ。財1を選ぶ確率は

\[ P(u_{i1} > u_{i0}) = P(\varepsilon_{i0} - \varepsilon_{i1} < v_{i1} - v_{i0}) \]

ここで、2つの独立なGumbel確率変数の差 \(\varepsilon_{i0} - \varepsilon_{i1}\) は、ロジスティック分布に従うという性質が使える（Gumbel分布の重要な性質の1つ）。ロジスティック分布の累積分布関数は \(\Lambda(x) = 1/(1+\exp(-x))\) なので、

\[ P(u_{i1} > u_{i0}) = \Lambda(v_{i1} - v_{i0}) = \frac{1}{1 + \exp(-(v_{i1}-v_{i0}))} = \frac{\exp(v_{i1})}{\exp(v_{i0}) + \exp(v_{i1})} \]

分母分子を \(\exp(v_{i0})\) で割れば、これは先ほどのlogit公式の \(J=1\) の場合そのものである。

\(J\)財版（max-stabilityを使う一般証明）

\(J+1\)個の選択肢がある一般の場合を考える。財 \(j\) を選ぶ確率は

\[ s_{ij} = P(u_{ij} \geq u_{ik}, \ \forall k) = P(\varepsilon_{ik} \leq v_{ij} - v_{ik} + \varepsilon_{ij}, \ \forall k \neq j) \]

\(\varepsilon_{ij} = \varepsilon\) の値で条件付けると、\(\varepsilon_{ik}\) が互いに独立なので、

\[ P(\varepsilon_{ik} \leq v_{ij} - v_{ik} + \varepsilon \ \forall k \mid \varepsilon_{ij}=\varepsilon) = \prod_{k \neq j} \exp\left(-\exp(-(v_{ij}-v_{ik}+\varepsilon))\right) \]

これを \(\varepsilon_{ij}\) の密度 \(f(\varepsilon) = \exp(-\varepsilon)\exp(-\exp(-\varepsilon))\) で積分すると（計算はやや込み入るが、指数関数の性質だけで閉じる）、

\[ s_{ij} = \int_{-\infty}^{\infty} \exp(-\varepsilon)\exp\left(-\exp(-\varepsilon)\sum_{k=0}^{J}\exp(v_{ik}-v_{ij})\right) d\varepsilon = \frac{\exp(v_{ij})}{\sum_{k=0}^{J}\exp(v_{ik})} \]

という閉形式が得られる。この積分が閉じるのが、まさにGumbel分布のmax-stability（独立なGumbelの最大値が再びGumbelになるという性質）の帰結である。他の分布（正規分布など）ではこの積分は一般に閉じず、数値積分やシミュレーションが必要になる。

Rシミュ1：Gumbel乱数で公式を確かめる

シミュレーションしてみよう。

Gumbel分布に従う乱数は、一様分布 runif() から簡単に作れる。\(U \sim \text{Uniform}(0,1)\) のとき、\(-\log(-\log(U))\) はGumbel分布に従う（Gumbel分布の逆関数法によるサンプリング）。

J <- 3          # ブランドA, B, C（+ outside option）
N <- 100000     # シミュレーションの消費者数

# 真の決定的効用（outsideを0に正規化）
delta_true <- c(0.8, 0.3, -0.2)   # ブランドA, B, C

# Gumbel乱数を引く: -log(-log(U))
draw_gumbel <- function(n) {
  u <- runif(n)
  -log(-log(u))
}

# N人 x (J+1)個の選択肢 の効用を作る
eps_mat <- matrix(draw_gumbel(N * (J + 1)), nrow = N, ncol = J + 1)
delta_full <- c(0, delta_true)  # outside=0を含めたJ+1本のベクトル

u_mat <- sweep(eps_mat, 2, delta_full, "+")
choice <- apply(u_mat, 1, which.max) - 1  # 0=outside, 1..J=ブランド

emp_freq <- sapply(0:J, function(j) mean(choice == j))

# logit理論式による選択確率
denom <- 1 + sum(exp(delta_true))
theory_prob <- c(1 / denom, exp(delta_true) / denom)

compare_df <- data.frame(
  alt = c("outside", "A", "B", "C"),
  empirical = round(emp_freq, 4),
  logit_formula = round(theory_prob, 4)
)
knitr::kable(compare_df, caption = "シミュレーションによる選択頻度 vs logit公式（N=100,000）")

シミュレーションによる選択頻度 vs logit公式（N=100,000）

alt	empirical	logit_formula
outside	0.1832	0.1854
A	0.4137	0.4126
B	0.2495	0.2502
C	0.1536	0.1518

plot_df <- compare_df %>%
  pivot_longer(cols = c(empirical, logit_formula), names_to = "type", values_to = "prob")

ggplot(plot_df, aes(x = alt, y = prob, fill = type)) +
  geom_col(position = position_dodge(width = 0.7), width = 0.6) +
  scale_fill_manual(values = c("empirical" = "#4C72B0", "logit_formula" = "#DD8452"),
                     labels = c("シミュレーション頻度", "logit理論確率")) +
  labs(x = "選択肢", y = "確率", fill = "",
       title = "Gumbel乱数からの選択 vs logit公式") +
  theme(legend.position = "top")

図 1: 経験選択頻度とlogit理論確率の比較（N=100,000）

見ての通り、2本の棒はほぼ完全に重なる。N=100,000まで増やせば、シミュレーションで得た選択頻度は理論値に収束していく。

これが、logit公式が単なる「魔法の数式」ではなく、Gumbel分布を仮定したRUMから機械的に導かれる帰結であることの、体感を伴った確認である。

これはどうでもいいことのようで、経済学的には結構大事である。Logit選択確率を使ういわゆるLogit modelは計量経済学や実証IO以外にも無数の科学分野で歴史上幾度となく発見されてきた。それをなぜ今更こんなにありがたがっているのか？それは、logit選択確率に「経済学的な解釈」がつけられたからである。error termがGumbel分布に従うときの効用最大化は側から見るとLogit 選択確率に従ったランダムな行動に見える、というある種のmicrofoudnationがここにはある。経済学者はこういうのが大好物である。天下り式に「ここは人々がlogitに従って選択しているとします」という話だったら何にも面白くない。「ここは人々のrandom utilityがGumbel distributionに従っていて、効用最大化で行動しています」となる、表面上は全く同じなのにとても嬉しくなる。そういう学問である。

Outside option

先ほどから当たり前のように使っている \(j=0\)（outside option）について、少し立ち止まって整理する。

離散選択モデルでは、効用のレベルには意味がなく、差にしか意味がない（\(u_{ij} - u_{ik}\) だけが選択確率に効く）。これはなぜかというと、 logitの式をoutside optionを基準に書き直すと、効用の差だけが残るからである。まだ \(v_{i0}=0\) と置かずに、分母・分子を \(\exp(v_{i0})\) で割ってみよう。

\[ \begin{aligned} s_{ij} &=\frac{\exp(v_{ij})}{\exp(v_{i0})+\sum_{k=1}^{J}\exp(v_{ik})}\\ &=\frac{\exp(v_{ij})/\exp(v_{i0})} {1+\sum_{k=1}^{J}\exp(v_{ik})/\exp(v_{i0})}\\ &=\frac{\exp(v_{ij}-v_{i0})} {1+\sum_{k=1}^{J}\exp(v_{ik}-v_{i0})}. \end{aligned} \]

選択確率そのものをoutsideの選択確率 \(s_{i0}\) で割っても、共通の分母が消えて、

\[ \frac{s_{ij}}{s_{i0}} =\frac{\exp(v_{ij})}{\exp(v_{i0})} =\exp(v_{ij}-v_{i0}) \]

となる。つまり、財 \(j\) が「買わない」よりどれだけ選ばれやすいかを決めるのは、\(v_{ij}\) の絶対的な大きさではなく、outsideとの差 \(v_{ij}-v_{i0}\) である。

全選択肢の効用に同じ定数を足しても選択は変わらないので、慣習としてoutside option（\(j=0\)）の決定的効用を \(v_{i0}=0\) と置く。これが位置（location）の正規化である。

もう1つ、効用全体を正の定数倍しても最大の選択肢は変わらない。誤差項を一般のスケール \(\sigma_\varepsilon\) のGumbel分布とすると、logit確率に現れるのは \(v_{ij}/\sigma_\varepsilon\) であり、データから識別されるのは係数そのものではなく「誤差のスケールに対する係数」である。

具体的には、\(\sigma_\varepsilon>0\) とし、標準Gumbel分布に従う独立な乱数を \(z_{ij}\) と書けば、

\[ \varepsilon_{ij}=\sigma_\varepsilon z_{ij}, \qquad u_{ij}=v_{ij}+\sigma_\varepsilon z_{ij}. \]

全選択肢の効用を同じ正の数 \(\sigma_\varepsilon\) で割っても、大小関係は変わらない。したがって、

\[ u_{ij}\geq u_{ik} \iff \frac{v_{ij}}{\sigma_\varepsilon}+z_{ij} \geq \frac{v_{ik}}{\sigma_\varepsilon}+z_{ik}. \]

右辺は、決定的効用が \(v_{ij}/\sigma_\varepsilon\)、誤差が標準Gumbelのモデルなので、先ほどのlogit公式をそのまま使えて、

\[ s_{ij} =\frac{\exp(v_{ij}/\sigma_\varepsilon)} {\sum_{k=0}^{J}\exp(v_{ik}/\sigma_\varepsilon)} \]

となる。例えば、任意の \(c>0\) について効用全体を \(c\) 倍すると、決定的効用は \(cv_{ij}\)、誤差のスケールは \(c\sigma_\varepsilon\) になるが、

\[ \frac{\exp\!\left(\frac{cv_{ij}}{c\sigma_\varepsilon}\right)} {\sum_{k=0}^{J}\exp\!\left(\frac{cv_{ik}}{c\sigma_\varepsilon}\right)} = \frac{\exp(v_{ij}/\sigma_\varepsilon)} {\sum_{k=0}^{J}\exp(v_{ik}/\sigma_\varepsilon)} =s_{ij}. \]

決定的効用と誤差のスケールを一緒に変えると、選択確率は全く同じである。今回の \(v_{ij}=x_j'\beta-\alpha p_j\) を代入すれば、

\[ \frac{v_{ij}}{\sigma_\varepsilon} =x_j'\left(\frac{\beta}{\sigma_\varepsilon}\right) -\left(\frac{\alpha}{\sigma_\varepsilon}\right)p_j. \]

このため、選択のデータから区別できるのは \(\beta/\sigma_\varepsilon\) と \(\alpha/\sigma_\varepsilon\) であり、\(\beta,\alpha,\sigma_\varepsilon\) をそれぞれ別々に決めることはできない。

そこで標準logitは \(\sigma_\varepsilon=1\) と置く。これが尺度（scale）の正規化である。

不変性	何を変えても選択が同じか	この講義の正規化
位置	全効用に同じ定数を足す	\(v_{i0}=0\)
尺度	全効用を同じ正の定数倍する	Gumbel誤差のスケールを1

したがって、異なるデータやモデル間で効用係数の絶対値を比べるときは、誤差スケールも違い得ることに注意が必要である。

ただし、実務上はここに悩ましい判断が入る。「買わない」という選択肢に、具体的に誰を含めるかという問題である。缶コーヒー棚の例で言えば、

• その棚の前を通った人のうち、何も買わずに通り過ぎた人（母数＝来店客数、あるいは棚の前を通った人数）
• 缶コーヒー以外の飲料（ペットボトルのコーヒーなど）を買った人も「缶コーヒーを買わない」に含めるか
• そもそも「市場規模」をどう定義するか（1日の来店客数？ 商圏人口？）

この市場規模（market size）の設定次第で、outside optionのシェアが変わり、ひいては価格弾力性の計算結果も変わる。この講義では単純化のため「棚の前で3ブランドのどれかを買うか、買わないかを決める人がN人いる」という設定で進めるが、実務ではこの市場規模の定義そのものが分析の質を左右する重要な論点であることは覚えておいてほしい。

選択確率を需要関数として読む

ここで、logit公式を「需要関数」として読み直そう。全員が同じ商品属性・価格に直面するとし、非価格部分の効用を \(a_j=x_j'\beta\)、推定する係数をまとめて \(\theta=(\beta,\alpha)\) と書く。このとき、商品 \(j\) の選択確率と予測需要量は

\[ s_j(\mathbf p,\mathbf x;\theta) =\frac{\exp(a_j-\alpha p_j)}{1+\sum_{k=1}^{J}\exp(a_k-\alpha p_k)}, \qquad q_j(\mathbf p,\mathbf x;\theta)=M s_j(\mathbf p,\mathbf x;\theta). \]

\(M\) は1回に最大1個を買う潜在的な消費者の数である。推定に使う標本人数 \(N\) と区別するため、予測対象の市場規模は \(M\) と書いた。右辺に自分の価格だけでなく競合の価格も入っているので、これは価格の組合せを各商品の需要量に変換する需要関数そのものである。共通の確率を持たない消費者の場合も、各人の購入確率を足した \(q_j=\sum_{i=1}^{M}s_{ij}\) が予測需要量になる。

手で1回計算してみよう。商品Aとoutsideだけの市場で、\(a_A=2\)、\(\alpha=0.02\)、\(M=1{,}000\) とする。Aが100円なら、

\[ v_A=2-0.02\times100=0,\qquad s_A=\frac{e^0}{1+e^0}=0.5,\qquad q_A=500\text{ 本}. \]

属性への評価と価格係数を固定して、価格だけを110円に変えると、

\[ v_A=2-0.02\times110=-0.2,\qquad s_A=\frac{e^{-0.2}}{1+e^{-0.2}}\simeq0.4502,\qquad q_A\simeq450.2\text{ 本}. \]

約450.2本は平均的に売れると予測する本数であって、1人が0.2本を買うという意味ではない。このように価格を順に変えて \(q_A\) を計算すれば、需要曲線が描ける。

2財それぞれの需要曲線を描く

今度は、商品A・商品B・outsideから1つを選ぶ市場に広げる。\(M=1{,}000\)、\(\alpha=0.02\) はそのままに、非価格部分の魅力を \(a_A=2\)、\(a_B=1.5\) としよう。logitから導かれる需要関数は、AにもBにも1本ずつある。

\[ \begin{aligned} q_A(p_A,p_B)&=1{,}000\frac{e^{2-0.02p_A}}{1+e^{2-0.02p_A}+e^{1.5-0.02p_B}},\\ q_B(p_A,p_B)&=1{,}000\frac{e^{1.5-0.02p_B}}{1+e^{2-0.02p_A}+e^{1.5-0.02p_B}}. \end{aligned} \]

まずは経済学でおなじみの、横軸に需要量、縦軸にその財自身の価格を置く図にしてみよう。左はAの需要曲線で、Bの価格を100円に固定してAの価格だけを動かす。右は逆に、Aの価格を100円に固定してBの価格だけを動かす。それぞれの青い実線が、その計算結果である。オレンジの破線は、もう一方の財の価格を150円に固定し直した場合を表す。

# A・B・outsideの3択。市場規模と属性への評価は固定する。
demand_two <- function(p_A, p_B) {
  exp_A <- exp(2 - 0.02 * p_A)
  exp_B <- exp(1.5 - 0.02 * p_B)
  denominator <- 1 + exp_A + exp_B
  data.frame(q_A = 1000 * exp_A / denominator,
             q_B = 1000 * exp_B / denominator,
             q_0 = 1000 / denominator)
}

demand_price_grid <- seq(0, 250, length.out = 251)
demand_curves <- do.call(rbind, lapply(c(100, 150), function(rival_price) {
  rbind(
    data.frame(product = "Aの需要曲線", own_price = demand_price_grid,
               quantity = demand_two(demand_price_grid, rival_price)$q_A,
               rival_price = factor(rival_price, levels = c(100, 150))),
    data.frame(product = "Bの需要曲線", own_price = demand_price_grid,
               quantity = demand_two(rival_price, demand_price_grid)$q_B,
               rival_price = factor(rival_price, levels = c(100, 150)))
  )
}))

ggplot(demand_curves, aes(x = quantity, y = own_price,
                         color = rival_price, linetype = rival_price)) +
  geom_hline(yintercept = 100, color = "gray80", linetype = "dotted") +
  geom_line(linewidth = 1) +
  geom_point(data = subset(demand_curves, own_price == 100), size = 2.3) +
  facet_wrap(~ product, nrow = 1) +
  scale_color_manual(values = c("100" = "#4C72B0", "150" = "#DD8452"),
                     name = "もう一方の財の価格", labels = c("100円", "150円")) +
  scale_linetype_manual(values = c("100" = "solid", "150" = "dashed"),
                       name = "もう一方の財の価格", labels = c("100円", "150円")) +
  scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, 250)) +
  scale_y_continuous(limits = c(0, 250), breaks = seq(0, 250, 50)) +
  labs(x = "その財の需要量（本）", y = "その財自身の価格（円）") +
  theme(legend.position = "bottom")

図 2: logitが導く2財の需要曲線。横軸は各財の需要量、縦軸は各財自身の価格。青い実線は競合価格100円、オレンジの破線は150円。点は自分の価格も100円のとき。

どちらも右下がりだが、直線ではない。価格が効用に入り、それを指数関数で変換しているので、価格帯によって数量の変わり方が異なる。また、需要量は0本と市場規模1,000本の間に収まる。outsideもあるため、\(q_A+q_B\) が必ず1,000本になるわけではない。

両方100円のとき、Aの需要は約384本、Bは約233本である。Aが先ほどの1財の例の500本より少ないのは、Bも選べるようになったためである。A自身は100円のまま、Bだけを150円にすると、Aの需要は約450本に増える。左の図では、青い点からオレンジの点へ右に移る。自分の価格を変えると1本の曲線に沿って動き、競合の価格を変えると需要曲線そのものが移動する。AもBも、需要が自分の価格だけで決まっているわけではないことが見えてくる。

本当は2つの価格に依存する：需要曲面とその断面

そこで、\(p_A\) と \(p_B\) の両方を動かしたくなる。\(q_A(p_A,p_B)\) は2つの数を受け取って需要量を返す関数だから、需要量を高さに取れば、3次元の曲面として描ける。\(q_B(p_A,p_B)\) についても、別の曲面が1枚ある。

下の図では、床の2軸をAの価格・Bの価格、高さを需要量とした。先ほどの2次元図とは縦軸の置き方を変えていることに注意しよう。Aの図の曲面は \(q_A\)、Bの図は \(q_B\) であり、どちらも同じ2つの価格に依存する。両図の床の軸は同じ向きにそろえてあり、同じ価格の組合せでAとBの需要量を比べられる。Aの曲面上の青い線は \(p_B=100\)、Bの曲面上の青い線は \(p_A=100\) に固定した断面である。これを取り出して価格を縦軸に置き直すと、先ほどの青い需要曲線になる。オレンジの線は、競合価格を150円に固定した断面に対応する。

ドラッグして回転 · 左右の視点は連動します

Aの需要量（本）

断面ではBの価格を固定

Bの需要量（本）

断面ではAの価格を固定

競合価格100円 競合価格150円

ノート静止画とRコード

draw_demand_surfaces <- function() {
  old_par <- par(mfrow = c(1, 2), mar = c(2.5, 1.5, 2.5, 1.5),
                 oma = c(1.5, 0, 0, 0),
                 family = if (Sys.info()[["sysname"]] == "Darwin") "Hiragino Sans" else "")
  on.exit(par(old_par))
  price_axis <- seq(0, 250, length.out = 41)
  quantity_A <- outer(price_axis, price_axis,
                      function(pa, pb) demand_two(pa, pb)$q_A)
  quantity_B <- outer(price_axis, price_axis,
                      function(pa, pb) demand_two(pa, pb)$q_B)

  for (product in c("A", "B")) {
    quantity <- if (product == "A") quantity_A else quantity_B
    projection <- persp(
      price_axis, price_axis, quantity,
      theta = 40, phi = 25, expand = 0.65,
      xlab = "Aの価格（円）", ylab = "Bの価格（円）",
      zlab = "",
      zlim = c(0, 1000), ticktype = "detailed", nticks = 3,
      col = "#E1E7EC", border = "#AEBCC7", shade = 0.25,
      cex.axis = 0.75, cex.lab = 0.9,
      main = paste0(product, "の需要量（高さ、単位：本）"))
    for (rival_price in c(100, 150)) {
      pa <- if (product == "A") price_axis else rep(rival_price, length(price_axis))
      pb <- if (product == "A") rep(rival_price, length(price_axis)) else price_axis
      q <- demand_two(pa, pb)[[paste0("q_", product)]]
      lines(trans3d(pa, pb, q, pmat = projection),
            col = if (rival_price == 100) "#4C72B0" else "#DD8452",
            lty = if (rival_price == 100) 1 else 2, lwd = 3)
    }
  }
  mtext("青い実線：競合価格100円の断面　／　オレンジの破線：競合価格150円の断面",
        side = 1, outer = TRUE, line = 0, cex = 0.85)
}
draw_demand_surfaces()

図 3: 2つの価格に依存する需要曲面。左の高さはAの需要量、右はBの需要量。曲面上の青い実線・オレンジの破線は、前の図と同じ競合価格100円・150円の断面を表す。

Aの曲面は、Aの価格を高くする方向では低くなり、Bの価格を高くする方向では高くなる。Bの曲面では、その役割が逆になる。続いて定義する自己価格弾力性は、曲面上で「自分の価格だけを動かす方向」の需要の反応を、交差価格弾力性は「もう一方の価格だけを動かす方向」の反応を、それぞれ変化率で測るものである。

需要の反応を弾力性で測る

需要曲面の図では、需要の反応を見る方向が2つあった。続けて、その「自分の価格を動かす方向」と「競合の価格を動かす方向」の反応を、変化率の比として測ってみよう。まず小さな数値例で計算し、その計算を微分による定義につなげる。

自分の価格を変える：自己価格弾力性

図と同じく、A・B・買わない、という3つの選択肢がある1,000人の市場を考える。ただし、ここでは計算を追いやすくするため、非価格部分の魅力を変え、Aの価格が100円のときのシェアをAが0.4、Bが0.3、買わない人が0.3とする。価格係数が \(\alpha=0.02\)（価格の単位は円）のlogitモデルで、Aだけを101円に値上げし、ほかの価格・属性と市場規模を固定すると、次のように予測される。

選択肢	値上げ前の期待人数	値上げ後の期待人数
Aを買う	400人	約395.210人
Bを買う	300人	約302.395人
買わない	300人	約302.395人

期待人数なので、小数になってよい。Aの効用だけが \(-0.02\times1=-0.02\) 変わるので、値上げ後のAのシェアは \(0.4e^{-0.02}/(0.4e^{-0.02}+0.3+0.3)\approx0.395210\) となる。

Aの価格の変化率は \((101-100)/100=1\%\)、Aの需要数量の変化率は \((395.2097-400)/400\approx-1.1976\%\) である。したがって、

\[ \frac{\text{Aの需要数量の変化率}}{\text{Aの価格の変化率}} \approx\frac{-1.1976\%}{1\%}=-1.1976 \]

となる。自分の価格を1%上げたとき、自分の需要が何%変わるかを測りたい、というのが自己価格弾力性（own-price elasticity）の発想である。ここで計算した値は、値上げ前を基準にした有限の変化率の比である。正式な弾力性は、価格変化を限りなく小さくしたときの比として、次に定義する。

競合の価格が変わる：交差価格弾力性

同じ値上げをBの側から見よう。B自身の価格は変えていないが、Aが高くなったことでBの需要は300人から約302.395人に増える。その変化率は \((302.3951-300)/300\approx0.7984\%\) なので、

\[ \frac{\text{Bの需要数量の変化率}}{\text{Aの価格の変化率}} \approx\frac{0.7984\%}{1\%}=0.7984 \]

となる。競合の価格を1%上げたとき、自分の需要が何%変わるかを測るのが交差価格弾力性（cross-price elasticity）の発想である。この例では正の値であり、AとBは代替財である。これも有限の変化率の比であり、弾力性の近似値として読む。

微分による定義とlogitの式

需要曲面で見た2つの方向の傾きは、Aについては \(\partial q_A/\partial p_A\) と \(\partial q_A/\partial p_B\) という偏微分に対応する。これらは「円あたり何本」という単位なので、先ほどの変化率の比にするには、価格と需要量で換算する必要がある。なお、最初の2次元図は価格が縦軸なので、図上の傾きは \(\partial p_A/\partial q_A\) の側である。需要の価格への反応と軸の向きを取り違えないようにしよう。

なぜ、傾きを「数量と価格の比」で割ると変化率の比になるのだろうか。価格を動かす前の数量と価格を \(q_j,p_k\)、そこからの変化を \(\Delta q_j,\Delta p_k\) と書こう。数量を%表示の変化率に直すには、変化した本数を元の本数で割って100を掛ける。価格についても同じことをすると、

\[ \begin{aligned} \frac{\text{需要数量の変化率（\%）}}{\text{価格の変化率（\%）}} &=\frac{100\,\Delta q_j/q_j}{100\,\Delta p_k/p_k}\\ &=\frac{\Delta q_j}{\Delta p_k}\frac{p_k}{q_j} =\frac{\Delta q_j/\Delta p_k}{q_j/p_k}. \end{aligned} \]

%表示のための100は、分子と分母で約分される。 残った式は、「数量の変化を価格の変化で割った傾き」を「元の数量と価格の比 \(q_j/p_k\)」で割ったものになっている。傾きと \(q_j/p_k\) はどちらも「本／円」という単位なので、その比には単位がない。弾力性は%変化どうしの比なので、例えば \((-1.2\%)/(1\%)=-1.2\) のように、比そのものには%を付けない。

数値例でも確認しよう。先ほどは1人1本なので、Aの需要は400本、Bは300本である。Aの価格100円の1%は1円、Aの需要400本の1%は4本である。1円の値上げによる約4.7903本の減少は、数量の1%に当たる4本の約1.1976倍である。つまり、価格が1%上がると需要が約1.1976%減る。同じ計算を式で書けば、

\[ \frac{\Delta q_A/\Delta p_A}{q_A/p_A} \approx\frac{-4.7903/1}{400/100} =\frac{-4.7903}{4} \approx-1.1976. \]

交差価格の場合も手順は同じである。Bの需要300本の1%は3本なので、Aの1円の値上げによるBの約2.3951本の増加は、Bにとって約0.7984%の増加になる。

\[ \frac{\Delta q_B/\Delta p_A}{q_B/p_A} \approx\frac{2.3951/1}{300/100} =\frac{2.3951}{3} \approx0.7984. \]

ここで数量の基準は反応を見るBの300本、価格の基準は値上げしたAの100円である。交差価格弾力性では、この組合せを取り違えないことが大切である。

ここまでの \(\Delta q_j/\Delta p_k\) は、2つの価格の間で測った平均的な傾きである。ほかの価格・属性と市場規模を固定し、価格変化を限りなく小さくすると、その点での傾き \(\partial q_j/\partial p_k\) になる。したがって、正の価格と需要数量のもとでは、

\[ \lim_{\Delta p_k\to0} \frac{\Delta q_j/q_j}{\Delta p_k/p_k} =\frac{\partial q_j}{\partial p_k}\frac{p_k}{q_j} =\frac{\partial q_j/\partial p_k}{q_j/p_k}. \]

この極限を弾力性と呼ぶ。動かす価格が自分の価格か、他の財の価格かによって、次の2つを区別する。

\[ \underbrace{\eta_{jj}=\frac{\partial q_j}{\partial p_j}\frac{p_j}{q_j}}_{\text{自己価格弾力性}}, \qquad \underbrace{\eta_{jk}=\frac{\partial q_j}{\partial p_k}\frac{p_k}{q_j}}_{\text{交差価格弾力性}\ (k\neq j)}. \]

\(\eta_{jk}\) の最初の添字 \(j\) は「需要が変わる財」、2番目の添字 \(k\) は「価格を変える財」である。小さな価格変化に対して、\(\Delta q_j/q_j\approx\eta_{jk}\Delta p_k/p_k\) と予測できる。市場規模 \(M\) を固定しているので、\(q_j=M s_j\) より、上の定義の \(q_j\) をシェア \(s_j\) に置き換えても同じ値になる。

logitでは、この定義にシェアの式を代入すると、

\[ \eta_{jj}=-\alpha p_j(1-s_j), \qquad \eta_{jk}=\alpha p_k s_k\quad(k\neq j) \]

が得られる。先ほどの値上げ前の状態では、

\[ \eta_{AA}=-0.02\times100\times(1-0.4)=-1.2, \qquad \eta_{BA}=0.02\times100\times0.4=0.8. \]

つまり、Aの1%の値上げに対して、Aの需要は約1.2%減り、Bの需要は約0.8%増えると近似できる。先ほどの正確なlogit予測（約1.1976%減、約0.7984%増）と少し違うのは、微分が値上げ前の1点での傾きを測るためである。

ノート導出：logit弾力性の式

outside optionの効用を0に基準化すると、\(s_j=\exp(v_j)/(1+\sum_{k=1}^J\exp(v_k))\) である。\(v_j=\ldots-\alpha p_j\) なので、\(\partial v_j/\partial p_j=-\alpha\) である。商の微分公式を使うと、

\[ \frac{\partial s_j}{\partial p_j} =-\alpha\cdot\frac{\exp(v_j)\left(1+\sum_{k=1}^J\exp(v_k)\right)-\exp(v_j)^2}{\left(1+\sum_{k=1}^J\exp(v_k)\right)^2} =-\alpha s_j(1-s_j). \]

弾力性の定義 \(\eta_{jj}=(\partial s_j/\partial p_j)(p_j/s_j)\) を使えば、

\[ \eta_{jj}=-\alpha s_j(1-s_j)\cdot\frac{p_j}{s_j} =-\alpha p_j(1-s_j). \]

交差弾力性も同様に、\(\partial s_j/\partial p_k=\alpha s_j s_k\)（\(k\neq j\)）から、

\[ \eta_{jk}=\alpha s_j s_k\cdot\frac{p_k}{s_j}=\alpha p_k s_k \]

が得られる。

最尤法（Maximum Likelihood Estimation）

ここまでは価格係数や属性への好みを与えて、需要曲線を描き、弾力性を計算してきた。では、その係数をデータからどう求めればよいのだろうか。観測された価格・属性のもとで人々が実際に何を選んだかを使い、その選択をよく説明する係数を探す方法として、最尤法（maximum likelihood estimation, MLE）を導入する。最尤法は計量経済学IIで初めて出てくる推定方法なので、じっくり丁寧に進める。

計量経済学Iでは、「推定量とはデータをパラメータの推定値に変換する関数である」という言い方をした。OLS推定量は、\((X,Y)\) のデータを受け取って \(\hat\beta = (X'X)^{-1}X'Y\) という関数でパラメータの推定値を返す関数だった。最尤法も同じで、データを受け取って、パラメータの推定値を返す、別の種類の関数である。ただし、その「関数の作り方」の発想がOLSとは違う。

3人の購入データから始める：Bernoulliと二項分布

まずは商品が1つだけで、「買う／買わない」の2択を考える。同じ価格を見た3人の記録が「買った、買わなかった、買った」だったとしよう。この価格で買う確率を、全員共通の未知の数 \(p\) と置く。ここだけの \(p\) は購入確率であり、商品の価格 \(p_j\) とは別の記号である。

試しに \(p=0.6\) という候補を考えると、各人の観測された行動に割り当てる確率は次の通りである。

消費者	観測された行動	その行動が起きる確率（\(p=0.6\)）
1人目	買った（\(y_1=1\)）	\(0.6\)
2人目	買わなかった（\(y_2=0\)）	\(1-0.6=0.4\)
3人目	買った（\(y_3=1\)）	\(0.6\)

3人の選択が独立だと仮定すると、この記録全体が得られる確率は掛け算で求められる。

\[ L(0.6)=0.6\times0.4\times0.6=0.144. \]

同じ記録に対して \(p=0.5\) なら \(L(0.5)=0.5\times0.5\times0.5=0.125\) である。したがって、この2候補の中では \(p=0.6\) の方が、観測された記録をより起こりやすいものとして説明する。データを固定し、候補のパラメータを変えながらこの確率を比較したものが尤度（likelihood）であり、最も大きくする候補を選ぶのが最尤法である。 尤度は「\(p\) 自体が正しい確率」ではない。

対数を取ると、積が和になる。自然対数を \(\ln\) と書けば、

\[ \ell(0.6)=\ln0.6+\ln0.4+\ln0.6 \simeq-0.511-0.916-0.511=-1.938. \]

\(p=0.5\) の対数尤度は \(3\ln0.5\simeq-2.079\) なので、対数を取った後も \(p=0.6\) の方が大きい。対数は単調増加だから最大にする候補は変わらず、たくさんの小さい確率を掛け続ける計算も避けられる。

ここまでの3人を \(N\) 人に増やそう。\(y_i=1\) なら確率 \(p\)、\(y_i=0\) なら確率 \(1-p\) を拾う操作は、\(p^{y_i}(1-p)^{1-y_i}\) と書ける。これが1人分のBernoulli分布の確率である。独立な \(N\) 人について掛け合わせ、対数を取ると、先ほどの計算をそのまま一般化できる。

\[ L(p)=\prod_{i=1}^{N}p^{y_i}(1-p)^{1-y_i},\qquad \ell(p)=\sum_{i=1}^{N}\left[y_i\ln p+(1-y_i)\ln(1-p)\right]. \]

では、誰が買ったかを区別せず、「3人中2人が買った」という人数だけをデータにしたらどうだろう。「買う・買う・買わない」「買う・買わない・買う」「買わない・買う・買う」の3通りを足すので、\(p=0.6\) での確率は \(3\times0.144=0.432\) になる。一般に \(N\) 人中 \(x=\sum_i y_i\) 人が買った人数は二項分布（binomial distribution）に従い、

\[ L_{\mathrm{count}}(p)=\binom{N}{x}p^x(1-p)^{N-x},\qquad \ell_{\mathrm{count}}(p)=\ln\binom{N}{x}+x\ln p+(N-x)\ln(1-p). \]

個々の記録の尤度と人数だけの尤度は異なるが、組合せの数 \(\binom{N}{x}\) は \(p\) に依存しないので、最尤推定値は同じである。コインの「表／裏」も、この「買う／買わない」と同じ計算になる。次の図では20回中14回表が出た記録について、定数項を除く対数尤度を描く。

n_coin <- 20
x_coin <- 14
p_grid <- seq(0.01, 0.99, by = 0.01)
loglik_coin <- x_coin * log(p_grid) + (n_coin - x_coin) * log(1 - p_grid)

coin_df <- data.frame(p = p_grid, loglik = loglik_coin)
p_hat_coin <- x_coin / n_coin  # 解析的な最尤推定量: 標本比率そのもの

ggplot(coin_df, aes(x = p, y = loglik)) +
  geom_line(color = "#4C72B0", linewidth = 1) +
  geom_vline(xintercept = p_hat_coin, linetype = "dashed", color = "#C44E52") +
  annotate("text", x = p_hat_coin + 0.05, y = min(loglik_coin) + 2,
           label = paste0("p_hat = ", p_hat_coin), color = "#C44E52") +
  labs(x = "p（表が出る確率の候補値）", y = "対数尤度",
       title = "対数尤度を最大にするpを探す")

図 4: コイン投げの対数尤度: n=20回中x=14回表が出た場合

図を見ると、対数尤度は \(p = x/n = 0.7\) でピークを迎える。実際、コイン投げの場合は微分して0とおけば解析的に \(\hat p = x/n\) が最尤推定量だと分かる（これは標本の表の比率そのものであり、直感とも一致する）。しかし、ロジットモデルではこの後で見るように解析解が求まらないため、数値的にピークを探す必要が出てくる。

logitの対数尤度

同じ掛け算を、商品A・商品B・買わない（outside）の3択でやってみよう。価格係数や商品属性への好みをまとめたパラメータの候補を \(\theta\) とする。ある候補をlogit公式に入れると、3人が直面する価格・属性のもとで次の確率が出たとしよう。各行の確率の合計は1である。

消費者	Aの確率	Bの確率	outsideの確率	実際に選んだもの	尤度に使う確率
1人目	0.6	0.2	0.2	A	0.6
2人目	0.2	0.5	0.3	B	0.5
3人目	0.3	0.2	0.5	outside	0.5

価格・属性を条件として3人の選択が独立なら、観測された「A、B、outside」という記録の尤度と対数尤度は、

\[ L(\theta)=0.6\times0.5\times0.5=0.15,\qquad \ell(\theta)=\ln0.6+\ln0.5+\ln0.5\simeq-1.897. \]

選ばれなかった商品の確率を全部掛けるわけではない。1人目がAを選ぶ確率はすでに0.6なので、その人について掛けるのは0.6だけである。別の候補 \(\theta'\) で、実際に選ばれた商品の確率が順に0.4、0.4、0.3になったなら、尤度は \(0.4\times0.4\times0.3=0.048\)、対数尤度は約 \(-3.037\) である。このデータには最初の候補の方がよく当てはまる。実際の推定では、各セルの確率をばらばらに選ぶのではなく、全員に共通する価格係数などを動かし、logit公式が返す確率をまとめて変える。

これを一般の人数・選択肢数で書こう。\(y_{ij}\) を「消費者 \(i\) が財 \(j\) を選んだら1、そうでなければ0」、\(s_{ij}(\theta)\) をモデルが返す選択確率とする。各人について \(y_{ij}=1\) は1つだけだから、\(\prod_j s_{ij}(\theta)^{y_{ij}}\) は、上の表の最後の列を拾う操作になる。

\[ L(\theta)=\prod_{i=1}^{N}\prod_{j=0}^{J}s_{ij}(\theta)^{y_{ij}}. \]

対数を取れば、掛け算が足し算に変わる。

\[ \ell(\theta) = \sum_{i=1}^{N} \sum_{j=0}^{J} y_{ij} \ln s_{ij}(\theta) \]

つまり、3人の例で行った「実際に選んだ確率のlogを足す」を、\(\sum\) で短く書いただけである。この \(\ell(\theta)\) を最大にする \(\hat\theta\) を、最尤推定量として採用する。

コイン投げとの橋を、2択の1人だけで確認しよう。財1を選ぶ確率を \(s_i(\theta)\)、実際の選択を \(y_i\in\{0,1\}\) とすると、その人の尤度は

\[ L_i(\theta)=s_i(\theta)^{y_i}\{1-s_i(\theta)\}^{1-y_i}. \]

\(y_i=1\) なら対数尤度への寄与は \(\log s_i\)、\(y_i=0\) なら \(\log(1-s_i)\) である。多項選択の \(\sum_j y_{ij}\log s_{ij}\) は、この「実際に選んだ確率のlogだけを拾う」を \(J+1\)択に拡張したものにすぎない。

二項分布のときと同様に、logitの対数尤度もグラフで眺められる。この後、実際にパラメータを推定してから、図 6 で価格係数を動かしたときの形を確認する。

線形指数のlogitモデルでは、対数尤度はパラメータについて大域的に凹（globally concave）である。したがって、見つかった局所最大は大域最大でもある。ただし、凹であることだけから「どの数値アルゴリズムも、どの初期値からでも必ず収束する」とは言えない。有限で一意な最大値には、説明変数行列の十分なランクと完全分離がないことが必要であり、平坦な方向・極端な初期値・数値誤差があれば最適化は失敗し得る。

Rシミュ2：自作MLEでlogitモデルを推定する

いよいよ、缶コーヒー棚の例に戻って、実際にデータをシミュレートし、自作の対数尤度を最大化してパラメータを推定してみる。

DGP（データ生成過程）

N=3,000人の消費者が、3つのブランド（A, B, C）+ outside optionから1つを選ぶ状況を考える。各ブランドには、価格（price）と甘さ（sweetness）という2つの属性があり、消費者は次の効用を持つとする。outsideを0に置く位置正規化に加え、ここではCのブランド固有切片を0とする追加の識別制約を置き、AとBだけに切片を置く。

\[ v_{ij} = \begin{cases} \text{const}_A + \beta_{\text{sweet}} \cdot \text{sweetness}_j - \alpha \cdot \text{price}_j & j = A \\ \text{const}_B + \beta_{\text{sweet}} \cdot \text{sweetness}_j - \alpha \cdot \text{price}_j & j = B \\ \beta_{\text{sweet}} \cdot \text{sweetness}_j - \alpha \cdot \text{price}_j & j = C \\ 0 & j = 0 \ (\text{outside}) \end{cases} \]

sweetnessはブランドごとに固定なので、A・B・Cすべての切片を入れると

\[ \text{sweetness}=3D_A+6D_B+5D_C \]

となり完全共線になる。現在の非価格部分のデザインは、行をoutside/A/B/C、列を \((\text{const}_A,\text{const}_B,\text{sweet})\) とすると

\[ \begin{pmatrix} 0&0&0\\ 1&0&3\\ 0&1&6\\ 0&0&5 \end{pmatrix} \]

で、列ランクは3である。Cの切片を0とするのはoutsideの位置正規化と同じではなく、「Cに残る観測不能なブランド魅力は0」という実質的な制約である。市場や商品をまたいでsweetnessが変化するデータがあれば、全ブランド切片と属性係数を同時に識別できる場合もある。 この3商品の例では、非価格部分は \(a_A=\text{const}_A+3\beta_{\text{sweet}}\)、\(a_B=\text{const}_B+6\beta_{\text{sweet}}\)、\(a_C=5\beta_{\text{sweet}}\) と書け、3つの \(a_j\) を自由に推定するのと同じだけの自由度がある。まずこの例で最尤法の実装を学び、後の6商品の例では、少数の共通の属性係数で商品間の効用を結びつける。

真のパラメータ値は次のように設定する。

N <- 3000

true_params <- list(
  const_A     = 1.2,
  const_B     = 0.6,
  beta_sweet  = 0.5,
  alpha_price = 0.02   # 価格は円単位。0.02 x 130円 ~ 2.6ほどの効用差になる程度のスケール
)

price_base <- c(130, 120, 110)   # ブランドA, B, C（円）
sweetness  <- c(3, 6, 5)          # 1-10スケール（甘さ）

価格係数 alpha_price を0.02としているのは恣意的な数字ではない。缶コーヒーの価格は100〜150円程度なので、alpha_price * price が1〜3程度の効用差になるようにスケールを合わせている。この手のスケール感覚は、シミュレーションを自分で書いて初めて身につくものなので、値を変えて挙動を見ることを宿題でやってもらう。

simulate_coffee_choice <- function(n, price_base, sweetness, params) {
  # 個人ごとにわずかな価格変動（プロモーションや店舗差）を与える
  price_noise <- matrix(rnorm(n * 3, mean = 0, sd = 8), nrow = n, ncol = 3)
  price <- sweep(price_noise, 2, price_base, "+")
  price <- pmin(pmax(price, 60), 200)

  delta <- cbind(
    params$const_A + params$beta_sweet * sweetness[1] - params$alpha_price * price[, 1],
    params$const_B + params$beta_sweet * sweetness[2] - params$alpha_price * price[, 2],
    0              + params$beta_sweet * sweetness[3] - params$alpha_price * price[, 3]
  )
  delta_full <- cbind(0, delta)  # outside = 0 の列を先頭に追加

  eps <- matrix(draw_gumbel(n * 4), nrow = n, ncol = 4)
  u <- delta_full + eps
  choice <- apply(u, 1, which.max) - 1  # 0=outside, 1=A, 2=B, 3=C

  list(choice = choice, price = price)
}

sim_data <- simulate_coffee_choice(N, price_base, sweetness, true_params)
choice_vec <- sim_data$choice
price_mat  <- sim_data$price

table(factor(choice_vec, levels = 0:3, labels = c("outside", "A", "B", "C")))

outside       A       B       C 
    439     490    1477     594 

share_df <- data.frame(
  alt = factor(c("outside", "A", "B", "C"), levels = c("outside", "A", "B", "C")),
  share = as.numeric(table(factor(choice_vec, levels = 0:3)) / N)
)

ggplot(share_df, aes(x = alt, y = share, fill = alt)) +
  geom_col(width = 0.6, show.legend = FALSE) +
  scale_fill_manual(values = c("outside" = "gray70", "A" = "#4C72B0", "B" = "#DD8452", "C" = "#55A868")) +
  labs(x = "選択肢", y = "シェア", title = "缶コーヒー棚：選択シェア")

図 5: 缶コーヒー3ブランド + outside の選択シェア（N=3,000）

おおよそ、outsideが14%、Bが最も甘く（sweetness=6）シェアも一番大きい（49%）、という直感的に妥当な分布になっている。

自作対数尤度とoptim()による最大化

対数尤度関数を自分で書く。\(\theta = (\text{const}_A, \text{const}_B, \beta_{\text{sweet}}, \alpha)\) の4パラメータをベクトルにまとめて受け取る関数として実装する。

build_design <- function(price, sweetness) {
  # 各消費者×各選択肢(outside,A,B,C)ごとの説明変数を配列で持つ
  # 列の並び: const_A, const_B, beta_sweet(の係数にかかる甘さ), alpha(の係数にかかる -price)
  n <- nrow(price)
  X <- array(0, dim = c(n, 4, 4))  # (消費者, 選択肢, パラメータ)
  # 選択肢1 = A
  X[, 2, 1] <- 1
  X[, 2, 3] <- sweetness[1]
  X[, 2, 4] <- -price[, 1]
  # 選択肢2 = B
  X[, 3, 2] <- 1
  X[, 3, 3] <- sweetness[2]
  X[, 3, 4] <- -price[, 2]
  # 選択肢3 = C（ブランド切片なし、Cが基準）
  X[, 4, 3] <- sweetness[3]
  X[, 4, 4] <- -price[, 3]
  X
}

neg_loglik_mnl <- function(theta, X, choice) {
  n <- dim(X)[1]
  J1 <- dim(X)[2]
  V <- matrix(0, nrow = n, ncol = J1)
  for (j in 1:J1) {
    V[, j] <- X[, j, ] %*% theta
  }
  V <- V - apply(V, 1, max)  # オーバーフロー対策（各行の最大値を引く。選択確率は変わらない）
  expV <- exp(V)
  denom <- rowSums(expV)
  chosen_idx <- cbind(1:n, choice + 1)  # choiceは0始まりなので+1
  logprob <- V[chosen_idx] - log(denom)
  -sum(logprob)
}

X_design <- build_design(price_mat, sweetness)

# optimは最小化するので、負の対数尤度を渡す
opt_result <- optim(
  par = c(0, 0, 0, 0),
  fn = neg_loglik_mnl,
  X = X_design,
  choice = choice_vec,
  method = "BFGS",
  hessian = TRUE
)

theta_hat <- opt_result$par
names(theta_hat) <- c("const_A", "const_B", "beta_sweet", "alpha_price")
theta_hat

    const_A     const_B  beta_sweet alpha_price 
 1.01343193  0.65368819  0.43280957  0.01699669 

true_theta <- c(true_params$const_A, true_params$const_B,
                 true_params$beta_sweet, true_params$alpha_price)

recovery_df <- data.frame(
  parameter = c("const_A", "const_B", "beta_sweet", "alpha_price"),
  true_value = true_theta,
  estimate = round(theta_hat, 4)
)
knitr::kable(recovery_df, caption = "真値とMLE推定値の比較（N=3,000）")

真値とMLE推定値の比較（N=3,000）

	parameter	true_value	estimate
const_A	const_A	1.20	1.0134
const_B	const_B	0.60	0.6537
beta_sweet	beta_sweet	0.50	0.4328
alpha_price	alpha_price	0.02	0.0170

推定値は、真値のごく近くに戻ってきている。これが「真値回収（parameter recovery）」の確認である。DGPを自分で書いてシミュレーションし、真のパラメータを与えてデータを作り、そのデータから同じパラメータを推定量が復元できるかどうかを確かめる——これは、推定量の実装が正しいかを検証する最も基本的で確実な方法である。実際のデータでは真値は分からないが、シミュレーションでは分かる。だからこそ、新しい推定量を実装したら、まずこの手順を踏むべきである。

logitの対数尤度も「山」として見る

二項分布の例では、確率 \(p\) を横軸にして対数尤度の山を描いた。今回のlogitでは4つの係数を同時に推定しているので、そのまま1枚の曲線にはできない。そこで、A・Bの切片と甘さの係数を推定値に固定し、価格係数 \(\alpha\) だけを動かした断面を描こう。3,000人の観測された価格・属性・選択はそのまま固定し、各候補の \(\alpha\) がその選択をどれだけよく説明するかを比べる。

alpha_mle <- unname(theta_hat["alpha_price"])
alpha_grid <- sort(unique(c(seq(0.005, 0.04, length.out = 201), alpha_mle)))

# 各候補でも同じデータを使う。ほかの3係数は推定値のままにする。
loglik_slice <- vapply(alpha_grid, function(alpha_try) {
  theta_try <- theta_hat
  theta_try["alpha_price"] <- alpha_try
  -neg_loglik_mnl(theta_try, X = X_design, choice = choice_vec)
}, numeric(1))
loglik_at_mle <- -neg_loglik_mnl(theta_hat, X = X_design, choice = choice_vec)
loglik_slice_df <- data.frame(alpha = alpha_grid, loglik = loglik_slice)

ggplot(loglik_slice_df, aes(x = alpha, y = loglik)) +
  geom_line(color = "#4C72B0", linewidth = 1) +
  geom_vline(xintercept = alpha_mle, linetype = "dashed", color = "#C44E52") +
  annotate("point", x = alpha_mle, y = loglik_at_mle,
           color = "#C44E52", size = 2.5) +
  annotate("text", x = alpha_mle + 0.001, hjust = 0,
           y = min(loglik_slice) + 0.1 * diff(range(loglik_slice)),
           label = paste0("alpha_hat = ", sprintf("%.4f", alpha_mle)),
           color = "#C44E52") +
  labs(x = "価格係数 α の候補値（価格は円単位）", y = "対数尤度",
       title = "対数尤度を最大にする価格係数を探す",
       subtitle = "ほかの3係数は推定値に固定")

図 6: logitの対数尤度の断面：価格係数だけを動かし、ほかの係数と観測データを固定

この断面でも、対数尤度は \(\hat\alpha=0.0170\) 付近で最大になる。二項分布では確率 \(p\) を直接動かしたが、logitでは係数 \(\alpha\) を動かすと各人の選択確率が変わり、その確率から対数尤度を計算する。観測データを固定して、対数尤度を最大にするパラメータを探すという考え方は同じである。ここで動かすのは商品の価格ではなく「価格にどれくらい反応するか」という係数である。また、各 \(\alpha\) でほかの係数を再推定するプロファイル尤度ではなく、推定値を通る固定した断面である。

ノート大域凹性について一言

logitの対数尤度は、\(\theta\) について大域的に凹である（ヘッセ行列が至る所で半負定値）。そのため、数値計算が収束して一階条件を満たした点は大域最大である。ただし、厳密な一意性にはデザイン行列の十分なランクと完全分離がないことが必要で、凹性はoptim()の収束自体を保証しない。上のコードで0を初期値にしたのは自然で安定しやすい選択だからであり、opt_result$convergence、勾配、別の初期値でも確認するのが実装上の作法である。

推論：HessianからのSEと信頼区間

推定値が得られただけでは、「本当にこの値なのか、それとも標本のブレでたまたまこの数字なのか」は分からない。計量経済学Iでは、OLS推定量の標準誤差（SE）から信頼区間や検定を作った。最尤推定量にも同じ発想が使える。

正しく特定された独立な選択の尤度のもとでは、\(\hat\theta\) は大標本で正規分布に近づき、分散共分散行列を負の対数尤度のHessian（観測情報行列）の逆行列で推定できる。optim()に hessian = TRUE を指定すると、最適化の最後に数値的なHessianを計算してくれる。ここでは負の対数尤度を最小化しているので、そのHessianをそのまま逆行列にする。

vcov_hat <- solve(opt_result$hessian)
se_hat <- sqrt(diag(vcov_hat))

ci_lower <- theta_hat - 1.96 * se_hat
ci_upper <- theta_hat + 1.96 * se_hat

inference_df <- data.frame(
  parameter = c("const_A", "const_B", "beta_sweet", "alpha_price"),
  true_value = true_theta,
  estimate = round(theta_hat, 4),
  se = round(se_hat, 5),
  ci_lower = round(ci_lower, 4),
  ci_upper = round(ci_upper, 4)
)
knitr::kable(inference_df, caption = "推定値・標準誤差・95%信頼区間（真値との比較）")

推定値・標準誤差・95%信頼区間（真値との比較）

	parameter	true_value	estimate	se	ci_lower	ci_upper
const_A	const_A	1.20	1.0134	0.20461	0.6124	1.4145
const_B	const_B	0.60	0.6537	0.06649	0.5234	0.7840
beta_sweet	beta_sweet	0.50	0.4328	0.06888	0.2978	0.5678
alpha_price	alpha_price	0.02	0.0170	0.00310	0.0109	0.0231

真値は、いずれのパラメータについても95%信頼区間の中に収まっている。推定値だけでなく、その不確実性の大きさまで含めて報告するという推論の作法は、計量経済学Iから変わらない。ただし、ここで計算した逆HessianはモデルベースのSEであり、ロバストSEそのものではない。尤度が誤特定されている場合はsandwich推定、同じ消費者の反復選択なら消費者単位のcluster-robust推定などが必要になる。このシミュレーションは正しいiid尤度から生成しているので、逆Hessianが適切である。

mlogitパッケージでの検算

自作の実装が正しいかを確かめるもう1つの方法は、信頼できる既存パッケージの結果と突き合わせることである。mlogitパッケージは離散選択モデル推定の標準的なRパッケージであり、ここでは検算のためだけに使う（本講義の主役はあくまで自作MLEである）。

library(mlogit)

# mlogitが要求するlong形式のデータフレームに変換する
mlogit_df <- data.frame(
  id = rep(1:N, each = 4),
  alt = rep(c("outside", "A", "B", "C"), times = N),
  # sapplyの各列が1人分(outside,A,B,C)。転置せず列方向に並べる
  chosen = as.vector(sapply(choice_vec, function(c) as.numeric(0:3 == c))),
  price = as.vector(t(cbind(0, price_mat))),
  sweet = rep(c(0, sweetness), times = N),
  asc_A = rep(c(0, 1, 0, 0), times = N),
  asc_B = rep(c(0, 0, 1, 0), times = N)
)

mlogit_data <- mlogit.data(mlogit_df, choice = "chosen", shape = "long",
                            alt.var = "alt", chid.var = "id")

mlogit_fit <- mlogit(chosen ~ asc_A + asc_B + sweet + price | 0, data = mlogit_data,
                      reflevel = "outside")
summary(mlogit_fit)

chid.var="id"は「この4行が同じ人の1つの選択集合」を指定する。chosenは各人について4行のうちちょうど1つがTRUEになっていることを必ず確認する。asc_Aとasc_Bを明示的に作り、| 0でパッケージ既定の選択肢固有切片を抑えているので、自作モデルと同じく「A・B切片あり、C切片なし」になる。priceをそのまま入れた係数は、自作コードの正のalpha_priceと符号が逆（\(-\hat\alpha\)）になる。もしchosen ~ price + sweet | 1とすると、A・B・Cの切片が自動で入り、ブランド固定のsweetnessと完全共線になるため、同じモデルの検算にはならない。ここでは記法の暗記より、2つの実装で同じデザイン行列を比較しているかを確認することが目的である。

Rシミュ3：標本分布とモンテカルロ

信頼区間が「95%」を名乗るからには、本当に95%の確率で真値を捕まえているのかを、モンテカルロ実験で確認しておきたい。同じDGPから何度もデータを生成し直し、そのたびに推定・信頼区間の計算を繰り返して、真値がどれくらいの頻度で信頼区間に入るかを数える。

n_mc <- 300
alpha_hat_mc <- numeric(n_mc)
covered_mc <- logical(n_mc)

for (m in 1:n_mc) {
  sim_m <- simulate_coffee_choice(N, price_base, sweetness, true_params)
  X_m <- build_design(sim_m$price, sweetness)

  opt_m <- optim(par = c(0, 0, 0, 0), fn = neg_loglik_mnl,
                 X = X_m, choice = sim_m$choice, method = "BFGS", hessian = TRUE)

  alpha_hat_mc[m] <- opt_m$par[4]  # alpha_priceの推定値だけ記録

  vcov_m <- tryCatch(solve(opt_m$hessian), error = function(e) NULL)
  if (!is.null(vcov_m)) {
    se_m <- sqrt(diag(vcov_m))[4]
    lo_m <- opt_m$par[4] - 1.96 * se_m
    hi_m <- opt_m$par[4] + 1.96 * se_m
    covered_mc[m] <- (lo_m <= true_params$alpha_price) && (true_params$alpha_price <= hi_m)
  } else {
    covered_mc[m] <- NA
  }
}

coverage_rate <- mean(covered_mc, na.rm = TRUE)
coverage_rate

[1] 0.9233333

mc_df <- data.frame(alpha_hat = alpha_hat_mc)
mc_mean <- mean(alpha_hat_mc)
mc_sd <- sd(alpha_hat_mc)

ggplot(mc_df, aes(x = alpha_hat)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "#4C72B0", alpha = 0.7) +
  stat_function(fun = dnorm, args = list(mean = mc_mean, sd = mc_sd),
                color = "#C44E52", linewidth = 1) +
  geom_vline(xintercept = true_params$alpha_price, linetype = "dashed", color = "black") +
  labs(x = expression(hat(alpha)), y = "density",
       title = paste0("alpha_priceの標本分布（モンテカルロ300回、被覆率=", round(coverage_rate, 3), "）"))

図 7: alpha_priceの標本分布（モンテカルロ300回）と正規近似

ヒストグラムは正規分布によく似た形になり、95%信頼区間の被覆率もおおむね0.95に近い値になる（乱数の具体的な引かれ方によって、0.94〜0.97程度の間で上下する）。この「モデルを書く→シミュレーションする→推定する→それを何度も繰り返して推定量の性質を確かめる」という一連の流れは、計量経済学IIを通して何度も使う型である。

Rシミュ4：商品の属性への係数を推定する

最尤法の手順が分かったので、今度は冒頭の特性アプローチをそのまま推定に使おう。苦味や甘さへの評価を推定できれば、「Aという商品がどれくらい人気か」に加えて、「Aの甘さを変えると、どれくらい売れるか」を考えられる。ここでは先ほどの3商品とは別のシミュレーションとして、6種類の缶コーヒーを考える。全商品に共通する苦味・甘さの係数を、消費者の選択から推定する。

属性から6商品の効用を作る

消費者はA〜Fの6商品とoutsideから1つを選ぶ。苦味と甘さは商品の属性を表すスコアであり、商品の組合せを次のようにする。Aの「苦味2・甘さ0」とBの「苦味0・甘さ2」は、冒頭の数値例と同じ組合せである。

products_attr <- data.frame(
  product = LETTERS[1:6],
  bitterness = c(2, 0, 0, 0, 2, 2),
  sweetness = c(0, 2, 0, 1, 1, 2),
  price_base = c(100, 100, 90, 95, 110, 120)
)
knitr::kable(products_attr,
             col.names = c("商品", "苦味", "甘さ", "基準価格（円）"))

商品	苦味	甘さ	基準価格（円）
A	2	0	100
B	0	2	100
C	0	0	90
D	0	1	95
E	2	1	110
F	2	2	120

全員に共通する属性の係数を使い、効用のうち観測した価格・属性で説明する部分を

\[ v_{ij}=\beta_0 +\beta_{\text{bitter}}\,\text{bitterness}_j +\beta_{\text{sweet}}\,\text{sweetness}_j -\alpha p_{ij}, \qquad j=A,\ldots,F \]

とする。\(p_{ij}\) は消費者 \(i\) に提示された商品 \(j\) の価格である。\(\beta_0\) は、どの商品にも共通する「outsideに対して商品を買うこと」の効用の部分であり、outsideの \(v_{i0}=0\) には加えない。たとえばAとBなら、

\[ \begin{aligned} v_{iA}&=\beta_0+2\beta_{\text{bitter}}-\alpha p_{iA},\\ v_{iB}&=\beta_0+2\beta_{\text{sweet}}-\alpha p_{iB}. \end{aligned} \]

同じ苦味が1単位増えれば、どの商品でも共通部分の効用は \(\beta_{\text{bitter}}\) だけ増える。この例では、6商品の非価格部分を、共通切片・苦味・甘さの3係数で説明する。商品ごとに独立したブランド切片は置かず、属性で説明しきれない商品固有の恒常的な魅力はないと仮定する。個人ごとの観測できない違いは、これまでと同じ独立なGumbel項 \(\varepsilon_{ij}\) によって表す。

真の係数を \(\beta_0=0.5\)、\(\beta_{\text{bitter}}=0.8\)、\(\beta_{\text{sweet}}=0.5\)、\(\alpha=0.02\) として、4,000人のデータを作る。価格には、各商品の基準価格の周りで個人ごとに独立な変動を与える。苦味・甘さの異なる商品間の選択に加え、同じ商品でも提示価格が異なると選択がどう変わるかが、係数を推定する手がかりになる。

set.seed(2026204)  # この6商品の例を独立に再現するためのseed
N_attr <- 4000
J_attr <- nrow(products_attr)
theta_attr_true <- c(beta0 = 0.5, beta_bitter = 0.8,
                     beta_sweet = 0.5, alpha_price = 0.02)

price_attr <- sweep(matrix(rnorm(N_attr * J_attr, sd = 8), nrow = N_attr),
                    2, products_attr$price_base, "+")
price_attr <- pmin(pmax(price_attr, 60), 200)

# outsideの行はすべて0。商品の行は (1, 苦味, 甘さ, -価格)。
build_attribute_design <- function(price, products) {
  n <- nrow(price)
  J <- nrow(products)
  X <- array(0, dim = c(n, J + 1, 4))
  for (j in seq_len(J)) {
    X[, j + 1, 1] <- 1
    X[, j + 1, 2] <- products$bitterness[j]
    X[, j + 1, 3] <- products$sweetness[j]
    X[, j + 1, 4] <- -price[, j]
  }
  X
}

X_attr <- build_attribute_design(price_attr, products_attr)
V_attr <- matrix(0, nrow = N_attr, ncol = J_attr + 1)
for (j in seq_len(J_attr + 1)) {
  V_attr[, j] <- X_attr[, j, ] %*% theta_attr_true
}
eps_attr <- matrix(draw_gumbel(N_attr * (J_attr + 1)), nrow = N_attr)
choice_attr <- max.col(V_attr + eps_attr, ties.method = "first") - 1

推定に使うのは、各人に提示された全商品の価格・属性と、実際に選んだ商品である。効用そのものやGumbel項は観測しない。データを作るために使った真の係数も、推定の際には伏せる。

選択データから共通の係数を推定する

係数の候補 \(\theta=(\beta_0,\beta_{\text{bitter}},\beta_{\text{sweet}},\alpha)\) を1組選ぶと、全商品の効用 \(v_{ij}(\theta)\) が同時に決まる。それをlogit公式に入れれば、

\[ s_{ij}(\theta) =\frac{\exp\{v_{ij}(\theta)\}} {1+\sum_{k=A}^{F}\exp\{v_{ik}(\theta)\}} \]

が各商品の選択確率となる。outsideの選択確率は、同じ分母のもとで分子を1としたものである。あとは、実際に選ばれた商品の対数確率を足した対数尤度を最大にすればよい。たとえば苦味の係数を大きくすると、Aだけでなく、同じく苦いEやFの予測効用も上がる。1つの属性係数の候補を、複数商品の選択データを使って評価しているのである。

前の推定と共通なのは、logit確率を計算して対数尤度を足す部分である。商品の効用を作る説明変数を \((1,\text{苦味},\text{甘さ},-\text{価格})\) に変えれば、先ほどの neg_loglik_mnl() をそのまま使える。

opt_attr <- optim(
  par = c(0, 0, 0, 0), fn = neg_loglik_mnl,
  X = X_attr, choice = choice_attr,
  method = "BFGS", hessian = TRUE,
  control = list(parscale = c(1, 1, 1, 0.01),
                 maxit = 1000, reltol = 1e-11)
)
stopifnot(opt_attr$convergence == 0)
theta_attr_hat <- setNames(opt_attr$par, names(theta_attr_true))
se_attr <- sqrt(diag(solve(opt_attr$hessian)))

knitr::kable(data.frame(
  coefficient = c("購入の共通切片", "苦味の係数", "甘さの係数", "価格係数 α"),
  truth = as.numeric(theta_attr_true),
  estimate = as.numeric(theta_attr_hat),
  se = se_attr,
  lower = as.numeric(theta_attr_hat - 1.96 * se_attr),
  upper = as.numeric(theta_attr_hat + 1.96 * se_attr)),
  digits = 4,
  col.names = c("係数", "真値", "推定値", "標準誤差", "95%下限", "95%上限"),
  caption = "6商品の選択データから推定した、全商品に共通する属性の係数")

6商品の選択データから推定した、全商品に共通する属性の係数

係数	真値	推定値	標準誤差	95%下限	95%上限
購入の共通切片	0.50	0.6433	0.2068	0.2379	1.0487
苦味の係数	0.80	0.8246	0.0286	0.7685	0.8806
甘さの係数	0.50	0.4940	0.0303	0.4347	0.5533
価格係数 α	0.02	0.0213	0.0024	0.0167	0.0259

たとえば甘さの係数は約 0.494 と推定された。価格とほかの属性を固定して甘さを1単位上げると、効用の共通部分がこの値だけ上がる、という意味である。これはA〜Fのすべてに共通する係数であり、消費者ごとの好みの係数 \(\beta_i\) を推定したものではない。また、選択確率がその値だけ上がるわけでもない。何本余分に売れるかを知るには、推定した効用をもう一度logit公式に入れる必要がある。

ノート商品ごとの切片もすべて入れたらどうなるか

このデータでは、苦味と甘さは各商品に固定されている。A〜Fのダミーをすべて入れると、苦味も甘さもそのダミーの線形結合として書けるので、属性の係数と商品切片を別々には推定できない。今回は商品固有の切片を置かず、6商品の非価格部分を3係数で説明する制約を課している。苦味・甘さの組合せが異なり、価格にも属性と独立な変動があるため、その制約のもとで各係数を推定できる。共通購入切片をoutsideとの差として推定するためには、買わなかった人もデータに含める。

推定した係数で、商品の属性変更を予測する

Aの甘さだけを0から1へ変え、価格は100円に保つとしよう。ほかの商品も表の基準価格のままで、属性・商品数・市場規模を固定する。推定したAの効用は

\[ \begin{aligned} \hat v_A^{\text{変更前}} &=\hat\beta_0+2\hat\beta_{\text{bitter}}-100\hat\alpha,\\ \hat v_A^{\text{変更後}} &=\hat v_A^{\text{変更前}}+\hat\beta_{\text{sweet}}. \end{aligned} \]

となる。Aの効用が上がると、Aの選択確率が増えると同時に、ほかの商品やoutsideの確率も分母を通じて変わる。推定したモデルにそのまま代入し、1,000人の市場での予測数量を比べよう。

predict_attribute_shares <- function(theta, price, products) {
  utility <- theta["beta0"] + theta["beta_bitter"] * products$bitterness +
    theta["beta_sweet"] * products$sweetness - theta["alpha_price"] * price
  utility <- c(0, utility)  # outside
  exp_utility <- exp(utility - max(utility))
  setNames(exp_utility / sum(exp_utility), c("outside", products$product))
}

shares_attr_before <- predict_attribute_shares(
  theta_attr_hat, products_attr$price_base, products_attr)
products_attr_after <- products_attr
products_attr_after$sweetness[products_attr_after$product == "A"] <- 1
shares_attr_after <- predict_attribute_shares(
  theta_attr_hat, products_attr$price_base, products_attr_after)

knitr::kable(data.frame(
  product = names(shares_attr_before),
  quantity_before = 1000 * as.numeric(shares_attr_before),
  quantity_after = 1000 * as.numeric(shares_attr_after),
  change = 1000 * as.numeric(shares_attr_after - shares_attr_before)),
  digits = 1,
  col.names = c("選択肢", "変更前の予測数量", "変更後の予測数量", "増減"),
  caption = "Aの甘さを0から1へ変えたときの予測（商品は本、outsideは人）")

Aの甘さを0から1へ変えたときの予測（商品は本、outsideは人）

選択肢	変更前の予測数量	変更後の予測数量	増減
outside	140.9	127.4	-13.5
A	165.8	245.6	79.9
B	85.6	77.4	-8.2
C	39.4	35.6	-3.8
D	58.1	52.5	-5.6
E	219.5	198.5	-21.0
F	290.8	262.9	-27.8

Aの予測数量は約 165.8 本から約 245.6 本へ増える。ここでは、Aの新しい人気度を別に推定し直したのではなく、全商品の選択から推定した甘さの係数を、Aの属性変更に当てはめた。これが、効用を商品の属性で表して推定する利点である。この予測は、甘さ以外の観測できない魅力や共通の選好係数も、変更によって変わらないという仮定のもとで行っている。次は、同じように価格を変えたときの需要を予測し、企業の利益最大化につなげる。

推定した需要関数をプライシングに使う

ここまでで、価格や属性への反応を表す係数をデータから推定できた。前半では係数を与えて需要曲線と弾力性を計算したが、ここからは推定値を使って「値上げしたら何本売れるか、客がどこへ流れるか」を予測し、価格設定につなげる。

推定値から弾力性を計算する

前半の弾力性の公式を振り返ると、自己価格弾力性は \(-\alpha p_j(1-s_j)\)、交差価格弾力性は \(\alpha p_k s_k\) である。ここでは最初に推定したA・B・Cの3ブランドのモデルに戻り、推定値 \(\hat\alpha\) と各価格での予測シェアを代入してみよう。

compute_shares <- function(theta, price, sweetness) {
  X_j <- build_design(price, sweetness)
  n <- dim(X_j)[1]; J1 <- dim(X_j)[2]
  V <- matrix(0, nrow = n, ncol = J1)
  for (j in 1:J1) V[, j] <- X_j[, j, ] %*% theta
  V <- V - apply(V, 1, max)
  expV <- exp(V)
  expV / rowSums(expV)
}

# 平均価格でシェアと弾力性行列を計算する（市場平均で評価）
avg_price_mat <- matrix(price_base, nrow = 1)
shares_avg <- compute_shares(theta_hat, avg_price_mat, sweetness)[1, ]
names(shares_avg) <- c("outside", "A", "B", "C")

alpha_est <- theta_hat["alpha_price"]
p_full <- c(NA, price_base)  # outsideには価格がない

elasticity_mat <- matrix(NA, nrow = 3, ncol = 3, dimnames = list(c("A","B","C"), c("A","B","C")))
for (j in 1:3) {
  for (k in 1:3) {
    if (j == k) {
      elasticity_mat[j, k] <- -alpha_est * p_full[j + 1] * (1 - shares_avg[j + 1])
    } else {
      elasticity_mat[j, k] <- alpha_est * p_full[k + 1] * shares_avg[k + 1]
    }
  }
}
knitr::kable(round(elasticity_mat, 4),
             caption = "弾力性行列（行=シェアが変化する財、列=価格が変化する財）")

弾力性行列（行=シェアが変化する財、列=価格が変化する財）

	A	B	C
A	-1.8500	1.0058	0.3687
B	0.3596	-1.0338	0.3687
C	0.3596	1.0058	-1.5009

この表を見ると、Bの価格を動かしたときの交差弾力性は \(\eta_{A,B}=\eta_{C,B}\) になっている。ここでは「率」と「量」を分けることが大切である。\(j\neq k\)について、

見る量	logitでの式	値上げ先\(k\)を固定したとき
相対変化（交差弾力性）	\(\partial\log s_j/\partial\log p_k=\alpha p_k s_k\)	流入先\(j\)によらず同じ
絶対変化（シェアポイント）	\(\partial s_j/\partial p_k=\alpha s_j s_k\)	元の\(s_j\)が大きい先ほど大きい

したがって、AとCの増加率は同じだが、流入するシェアポイントまで同じとは限らない。絶対的な流出は各行き先の元のシェアに比例配分される。AがBに似ていても、その「近さ」による上乗せはなく、Cとの違いは元の人気度だけである。これがIIA（Independence of Irrelevant Alternatives）が課す代替パターンの制約である。

需要予測をプライシングにつなげる

需要関数が分かると、候補となる価格ごとに販売数量を予測できる。売り手が選びたいのは、その予測を使って利益を最大にする価格である。そこで、まずAだけを販売する1社の価格設定を考え、次にAとBを別々の企業が販売する価格競争へ進もう。どちらも市場規模と商品の属性は固定し、売り手は予測需要量をすべて供給できるとする。

1社がAだけを販売する：利益最大化と自己価格弾力性

消費者がAを買うか、何も買わないかを選ぶ市場を考える。売り手は需要関数 \(q_A(p_A)\) を把握しており、Aを1本追加で供給する費用（限界費用）は一定の \(c_A\) 円である。1本売ると \(p_A-c_A\) 円が残るので、価格で変わらない固定費を省いた利益は \((p_A-c_A)q_A(p_A)\) となる。売り手の問題は

\[ \max_{p_A\geq 0}\ \pi_A(p_A) =\max_{p_A\geq 0}\ (p_A-c_A)q_A(p_A) \]

である。価格を上げると1本当たりの利幅は増えるが、売れる本数は減る。この2つの効果がちょうど釣り合う価格を探したい。需要が微分可能で、最適価格が正の内点にあるなら、利益の微分が0になるという一階条件は

\[ \frac{d\pi_A}{dp_A} =\underbrace{q_A(p_A)}_{\text{単価上昇による利益の増加}} +\underbrace{(p_A-c_A)\frac{dq_A(p_A)}{dp_A}}_{\text{数量減少による利益の減少}} =0 \]

となる。価格をわずかに、たとえば1円引き上げると、もともと売れていた各1本から約1円ずつ余分に得られる。一方、買わなくなる人の分だけ、1本当たり \(p_A-c_A\) 円の利益を失う。

ここで自己価格弾力性 \(\eta_{AA}(p_A)=\{dq_A(p_A)/dp_A\}\,p_A/q_A(p_A)\) を使うと、この条件を読みやすく書き直せる。上式を正の需要量 \(q_A(p_A)\) で割り、微分の項を弾力性に置き換えると、

\[ 1+\frac{p_A-c_A}{p_A}\eta_{AA}(p_A)=0. \]

したがって、最適価格 \(p_A^*\) は

\[ \boxed{\frac{p_A^*-c_A}{p_A^*} =-\frac{1}{\eta_{AA}(p_A^*)} =\frac{1}{|\eta_{AA}(p_A^*)|}} \]

を満たす。価格に占める利幅の割合が、自己価格弾力性の絶対値の逆数に等しくなるという関係である。値上げで数量が大きく減る需要に直面するほど、最適価格で維持できる利幅の割合は小さい。ただし、弾力性自体も価格とともに変わる。今の価格で測った弾力性を固定して値付けするのではなく、需要と弾力性を計算し直しながら、この条件を満たす価格を求める。

前半のAとoutsideだけの数値例に戻ろう。\(M=1{,}000\)、\(a_A=2\)、\(\alpha=0.02\) なので、

\[ s_A(p_A)=\frac{e^{2-0.02p_A}}{1+e^{2-0.02p_A}}, \qquad q_A(p_A)=1{,}000s_A(p_A). \]

ここに \(c_A=50\) 円という費用の条件を加える。logitの自己価格弾力性 \(\eta_{AA}=-\alpha p_A(1-s_A)\) を一階条件に代入すると、解くべき式は

\[ p_A^*=50+\frac{1}{0.02\{1-s_A(p_A^*)\}} \]

となる。右辺のシェアにも \(p_A^*\) が入っているので、Rで数値的に解いてみよう。

pricing_alpha <- 0.02
pricing_cost <- 50
share_single <- function(p_A) plogis(2 - pricing_alpha * p_A)

# 一階条件 1 - alpha * (p_A - c_A) * (1 - s_A) = 0 を解く。
price_single <- uniroot(function(p_A) {
  1 - pricing_alpha * (p_A - pricing_cost) * (1 - share_single(p_A))
}, interval = c(pricing_cost, 500), tol = 1e-10)$root
quantity_single <- 1000 * share_single(price_single)
profit_single <- (price_single - pricing_cost) * quantity_single
eta_single <- -pricing_alpha * price_single * (1 - share_single(price_single))

答えは約 128.36 円である。このとき約 361.9 本売れ、利益は約 28,357 円になる。弾力性は約 -1.638、価格に占める利幅は約 0.610 であり、上の逆数の関係も確認できる。

収入（売上金額）だけを最大にする場合は、目的関数が \(R_A(p_A)=p_Aq_A(p_A)\) になる。これは費用を0として同じ計算をすることに等しく、一階条件は \(\eta_{AA}=-1\) である。この例では収入最大化価格は100円で、500本が売れる。そのときの利益は \((100-50)\times500=25{,}000\) 円なので、利益最大化価格での利益より小さい。下の図では、収入と利益の頂点が異なる位置に現れる。

pricing_df <- data.frame(price = seq(50, 250, length.out = 401))
pricing_df$share <- share_single(pricing_df$price)
pricing_df$revenue <- pricing_df$price * 1000 * pricing_df$share
pricing_df$profit <- (pricing_df$price - pricing_cost) * 1000 * pricing_df$share
pricing_money <- data.frame(
  price = rep(pricing_df$price, 2),
  amount = c(pricing_df$revenue, pricing_df$profit),
  measure = factor(rep(c("収入", "利益"), each = nrow(pricing_df)),
                   levels = c("収入", "利益")))
pricing_optima <- data.frame(
  price = c(100, price_single),
  measure = factor(c("収入", "利益"), levels = c("収入", "利益")))

p1 <- ggplot(pricing_df, aes(x = price, y = share)) +
  geom_line(color = "#4C72B0", linewidth = 1) +
  labs(x = "Aの価格（円）", y = "予測シェア", title = "価格を上げると需要が減る")
p2 <- ggplot(pricing_money, aes(x = price, y = amount / 1000, color = measure)) +
  geom_line(linewidth = 1) +
  geom_vline(data = pricing_optima, aes(xintercept = price, color = measure),
             linetype = "dashed", show.legend = FALSE) +
  scale_color_manual(values = c("収入" = "#4C72B0", "利益" = "#DD8452"), name = NULL) +
  labs(x = "Aの価格（円）", y = "金額（千円）", title = "収入と利益では最適価格が違う") +
  theme(legend.position = "bottom")
p1 + p2

図 8: Aとoutsideだけの市場。限界費用は50円。左は予測シェア、右は収入と利益。縦線は各曲線を最大にする価格を示す。

AとBを別々の企業が販売する：最適反応と価格の均衡

次に、Aを売る企業とBを売る企業が1社ずついる市場を考える。各企業は自分の商品の価格を同時に選び、消費者はA・B・outsideから選ぶ。前半で描いた2つの需要曲面 \(q_A(p_A,p_B)\)、\(q_B(p_A,p_B)\) が、ここで必要になる。Aの販売数量はAの価格だけでなく、Bの価格にも依存するからである。

各商品の一定の限界費用を \(c_A\)、\(c_B\) とする。Aの企業は、Bの価格を所与として自分の利益を最大にする。同様にBの企業は、Aの価格を所与として自分の利益を最大にする。各社の問題は

\[ \begin{aligned} \text{Aの企業：}\quad &\max_{p_A\geq0}\ (p_A-c_A)q_A(p_A,p_B) &&\text{（$p_B$を固定）},\\ \text{Bの企業：}\quad &\max_{p_B\geq0}\ (p_B-c_B)q_B(p_A,p_B) &&\text{（$p_A$を固定）}. \end{aligned} \]

相手の価格に対して、自分の利益を最大にする価格を最適反応と呼ぶ。そして、互いの価格が相手に対する最適反応になっている価格の組 \((p_A^*,p_B^*)\) が、価格競争のナッシュ均衡である。均衡では、相手が価格を変えないかぎり、どちらの企業も自分だけ価格を変更して利益を増やせない。式では

\[ \begin{aligned} p_A^*&\in\operatorname*{arg\,max}_{p_A\geq0}(p_A-c_A)q_A(p_A,p_B^*),\\ p_B^*&\in\operatorname*{arg\,max}_{p_B\geq0}(p_B-c_B)q_B(p_A^*,p_B) \end{aligned} \]

と書ける。\(\operatorname*{arg\,max}\) は、最大の利益額そのものではなく、その利益を実現する価格を表す。

均衡価格を計算するにも、弾力性が役に立つ。各社が相手価格を固定して自分の価格で利益を微分すると、正の内点での一階条件は

\[ \begin{aligned} q_A+(p_A-c_A)\frac{\partial q_A}{\partial p_A}&=0,\\ q_B+(p_B-c_B)\frac{\partial q_B}{\partial p_B}&=0. \end{aligned} \]

1社の場合と同じ変形により、均衡では次の2つが同時に成り立つ。

\[ \frac{p_A^*-c_A}{p_A^*} =-\frac{1}{\eta_{AA}(p_A^*,p_B^*)}, \qquad \frac{p_B^*-c_B}{p_B^*} =-\frac{1}{\eta_{BB}(p_A^*,p_B^*)}. \]

どちらの企業も1商品だけを販売するので、自分の値上げでBやAに移った顧客から利益を得ることはない。そのため、それぞれの条件に直接入るのは自己価格弾力性である。一方、交差価格弾力性は、相手が価格を変えたときに自社の需要がどう変わるかを表す。その需要の変化によって自社の最適反応も変わるため、2社の価格を求めるには、それぞれの需要が両方の価格にどう依存するかを把握する必要がある。

数値例には、前半の需要曲面と同じ \(M=1{,}000\)、\(a_A=2\)、\(a_B=1.5\)、\(\alpha=0.02\) を使う。費用は両社とも1本50円とする。このとき、シェアと均衡の一階条件は

\[ \begin{aligned} s_A(p_A,p_B)&=\frac{e^{2-0.02p_A}}{1+e^{2-0.02p_A}+e^{1.5-0.02p_B}},\\ s_B(p_A,p_B)&=\frac{e^{1.5-0.02p_B}}{1+e^{2-0.02p_A}+e^{1.5-0.02p_B}},\\[6pt] p_A^*&=50+\frac{1}{0.02\{1-s_A(p_A^*,p_B^*)\}},\\ p_B^*&=50+\frac{1}{0.02\{1-s_B(p_A^*,p_B^*)\}}. \end{aligned} \]

今度は未知の価格が2つあり、右辺のシェアも両方の価格に依存する。この2つの条件を連立して解くと、次の均衡が得られる。

# 相手の価格を固定し、各社の一階条件を解く。
best_response_A <- function(p_B) {
  uniroot(function(p_A) {
    s_A <- demand_two(p_A, p_B)$q_A / 1000
    1 - pricing_alpha * (p_A - pricing_cost) * (1 - s_A)
  }, interval = c(pricing_cost, 500), tol = 1e-10)$root
}
best_response_B <- function(p_A) {
  uniroot(function(p_B) {
    s_B <- demand_two(p_A, p_B)$q_B / 1000
    1 - pricing_alpha * (p_B - pricing_cost) * (1 - s_B)
  }, interval = c(pricing_cost, 500), tol = 1e-10)$root
}

# 互いの最適反応を順に計算し、両価格が変わらなくなるまで繰り返す。
price_nash <- c(A = 100, B = 100)
for (iteration in 1:100) {
  previous_prices <- price_nash
  price_nash["A"] <- best_response_A(price_nash["B"])
  price_nash["B"] <- best_response_B(price_nash["A"])
  if (max(abs(price_nash - previous_prices)) < 1e-9) break
}
stopifnot(abs(price_nash["A"] - best_response_A(price_nash["B"])) < 1e-7,
          abs(price_nash["B"] - best_response_B(price_nash["A"])) < 1e-7)
quantity_nash <- unlist(demand_two(price_nash["A"], price_nash["B"])[c("q_A", "q_B")],
                        use.names = FALSE)
eta_nash <- -pricing_alpha * price_nash * (1 - quantity_nash / 1000)
knitr::kable(data.frame(
  商品 = c("A", "B"),
  均衡価格_円 = as.numeric(price_nash),
  予測需要量_本 = quantity_nash,
  自己価格弾力性 = as.numeric(eta_nash),
  利幅の割合 = as.numeric((price_nash - pricing_cost) / price_nash)),
  digits = c(0, 2, 1, 3, 3),
  col.names = c("商品", "均衡価格（円）", "予測需要量（本）", "自己価格弾力性", "利幅／価格"),
  caption = "AとBを別々の企業が販売する場合のナッシュ均衡（限界費用は両社とも50円）")

AとBを別々の企業が販売する場合のナッシュ均衡（限界費用は両社とも50円）

商品	均衡価格（円）	予測需要量（本）	自己価格弾力性	利幅／価格
A	122.05	306.1	-1.694	0.590
B	113.96	218.2	-1.782	0.561

Bが約 113.96 円なら、Aにとっては約 122.05 円が最も利益の高い価格である。同時に、Aがその価格なら、Bにとっては約 113.96 円が最も利益の高い価格になっている。Aだけを販売する市場では約 128.36 円だった最適価格が、この例ではBとの競争によって低くなる。

ノート一階条件を解けば、本当に利益最大になるのか

一般の需要関数では、微分が0になることは内点での必要条件であり、それだけで最大とは限らない。今回のlogit需要と一定の限界費用では、Aだけの市場でも、Bの価格を固定した市場でも、利益の微分は

\[ \frac{\partial\pi_A}{\partial p_A} =q_A\bigl\{1-\alpha(p_A-c_A)(1-s_A)\bigr\} \]

と書ける。\(p_A<c_A\) では利益が負なので、\(p_A\geq c_A\) を調べればよい。この範囲で価格を上げると、\(p_A-c_A\) は増え、\(s_A\) は下がるので \(1-s_A\) も増える。そのため、中括弧は \(p_A=c_A\) での1から単調に下がり、やがて負になる。利益は一階条件の解まで増え、その後は減るので、解は利益を最大にする唯一の価格である。Bについても同じ理由が成り立ち、両社の一階条件を同時に満たす上の価格は、互いに最適反応となっている。

この2つの例では、計算の流れが見えるように前半の係数を使った。推定した需要を使う場合も、\(\hat\alpha\) や属性の推定係数から各価格での需要量と弾力性を計算し、同じ利益最大化問題を解けばよい。その際には、需要の推定に加えて限界費用の情報が必要になる。競争相手も価格を変える場面では、さらに「各社が同時に自社の価格を選ぶ」といった競争の仕方を定めることで、予測需要から均衡価格まで計算できる。

実証研究コーナー

ヒント実証研究コーナー：McFaddenによるBART需要予測の事後検証

Daniel McFadden はmultinomial logitモデルの理論的基礎を築き、2000年にノーベル経済学賞を受賞した。McFadden (2001) の回顧によれば、研究チームはBART開業前の1972年に通勤行動を調査してMNLを推定し、開業後の1975年に再調査できた631人について、実現した各交通手段の属性を使って予測を評価した。BARTの予測シェアは6.3%、実績は6.2%だった（予測の標準誤差は2.5パーセントポイント）。一方で、徒歩アクセスを過大評価するなど系統的な誤差もあり、McFadden自身も精密な一致には幸運があったと注意している。したがってこれは「常に小数点まで当たる魔法」ではなく、個票RUMが当時の集計型公式予測（BARTシェア15%）より政策に敏感な予測を作れた、という検証事例として読むのが正確である。

ヒント実証研究コーナー：Guadagni and Little (1983)

Guadagni and Little (1983) “A Logit Model of Brand Choice Calibrated on Scanner Data” は、マーケティングサイエンスにおけるブランド選択モデルの古典的研究である。この論文が画期的だったのは、スーパーマーケットの実際のスキャナーデータを使ってコーヒーのブランド選択を推定し、価格・特売（feature）・陳列（display）に加えて、過去の購買を指数加重したブランドロイヤルティ変数を組み込んだ点にある。この変数は予測には強力だが、その係数を「過去に買ったから今も買う」という真のstate dependenceの因果効果と直ちに解釈してはいけない。もともと特定ブランドを好む持続的な未観測異質性も購買履歴に現れるため、両者が混ざり得る。この区別は第14回で改めて扱う。

ヒント実証研究コーナー：Train (2009) — 離散選択の標準教科書

Kenneth Train の教科書 Discrete Choice Methods with Simulation（2009年、第2版）は、離散選択モデルを学ぶ者にとって避けて通れない標準教科書である。この講義のReadingにも指定した第2〜3章は、まさに今日扱ったRUMとlogitモデルの導出を、非常に丁寧に扱っている。Trainのユニークな点は、書名の通り「シミュレーション」に重心を置いていることで、mixed logit（第4回で扱うrandom coefficient logitの発展形）のように解析的に選択確率が求まらないモデルについても、モンテカルロシミュレーションでどう推定するかを実践的に解説している。この講義で「まずシミュレーションで公式を確かめる」というスタイルを取っているのも、Trainの精神を受け継いだものである。オンラインでPDFが無料公開されていることでも知られており、本格的に離散選択モデルを実装する立場になったら、必ず手元に置いておきたい一冊である。

ヒント実証研究コーナー：動画配信のプラン選択

サブスクリプション型サービス（動画配信、音楽配信など）のプラン設計は、離散選択モデルの応用先として近年ますます重要度を増している。複数の価格帯・画質・同時視聴数などを組み合わせたプラン群から消費者がどれを選ぶかは、まさに今日学んだMNLの枠組みでそのまま扱える問題である。本講義のcoding課題（assignment2.qmd）でも、この設定を題材にする。実務では、新しいプランを追加したときに既存プランからどれだけ客が流れるか（カニバリゼーション）を事前に予測する必要があり、そこでもIIAの性質（今日の弾力性表で示唆した「交差弾力性が財に依らない」という制約）が問題になることがある。この点は次回の講義でさらに掘り下げる。

ビジネスの現場で

ノートビジネスの現場で：サブスクのプラン選択

Basic・Standard・Premium・未加入という4つの選択肢からの選択は、缶コーヒー棚の例とまったく同じ構造を持つ離散選択の問題である。月額料金、画質、同時視聴可能数といった属性が \(v_{ij}\) に入り、未加入が outside option になる。この枠組みが使えれば、「Standardプランを値上げしたら、どれだけの客がBasicにダウングレードし、どれだけが解約（outside）するか」という、サブスクビジネスで最も気になる問いに答える土台ができる。

ノートビジネスの現場で：POSデータから棚割り・値付けへ

コンビニやスーパーのPOSデータには、「いつ・どの店で・どの商品が・いくらで売れたか」という取引記録が大量に蓄積されている。単品ごとの売上を時系列で見るだけでは、値下げしたときにどの競合商品からシェアを奪ったのかは分からない。今日学んだように、同じ棚に並ぶ商品群を「選択肢の集合」として捉え、離散選択モデルを当てはめることで、価格弾力性や商品間の代替関係を推定でき、それが棚割り最適化や販促プランの設計に直結する。

重要ここが核心

• 製品を特徴の束として表すのが特性アプローチである。今回は属性と価格の係数を全員共通とし、観測できない個人ごとの違いを \(\varepsilon_{ij}\) で表す。個人の選択を集計すると市場の需要になる。
• 消費者の選択は、決定的効用 \(v_{ij}\) とランダム項 \(\varepsilon_{ij}\) からなるRUM \(u_{ij}=v_{ij}+\varepsilon_{ij}\) として書ける。
• \(\varepsilon_{ij}\) に独立なGumbel分布を仮定すると、選択確率がlogit公式になる。価格が効用に入るので、価格を変えると確率が変わり、その確率を人数分集計した \(q_j\) が需要関数になる。
• 最尤法は、「観測されたデータが最も得られやすいパラメータ」を選ぶ推定方法であり、logitの対数尤度 \(\sum_i\sum_j y_{ij}\ln s_{ij}\) を最大化する。効用を属性で表せば、共通の属性係数を推定し、商品の属性を変えたときの需要も予測できる。
• Hessianの逆行列から標準誤差・信頼区間が計算でき、モンテカルロ実験でその被覆率を検証できる。
• 自己価格弾力性は \(-\alpha p_j(1-s_j)\)、交差価格弾力性は \(\alpha p_k s_k\) という閉形式で書け、プライシングの意思決定に直結する。ただし交差弾力性が財の類似性を反映しない（IIA）という制約があり、これが次回のテーマになる。

宿題

• ブラウザ実験室（playground2.html）で、今回のシミュレーションの数値を自分でいじって遊んでみよう。インストール不要、ブラウザだけでRが動く（初回ロードのみ30秒〜1分かかる）。
• 今日のoptim()のコードで、method = "BFGS"をmethod = "Nelder-Mead"やmethod = "CG"に変えて実行し、推定値がどう変わるか（あるいは変わらないか）を確認してみよう。
• ChatGPTなどの生成AIに「なぜ離散選択モデルではGumbel分布（Type I極値分布）を使うのか、正規分布ではダメなのか」と聞いてみて、今日の講義の説明とどう違うか比べてみよう。
• coding課題は assignment2.qmd を参照。動画配信サービスのプラン選択を題材に、今日と同じ手順（DGP→自作MLE→推論→反実仮想）を自分の手で実装してもらう。

次回予告

次回（第3回）は、今日作ったMNLモデルの限界に切り込む。今日、弾力性行列の表で示唆した「交差弾力性が財の類似性を無視する」という性質は、IIA（Independence of Irrelevant Alternatives、無関係な選択肢からの独立性）と呼ばれる、logitモデルに内在する仮定の帰結である。この制約を緩和するnested logitモデルを導入し、さらに、価格が観測されない需要ショックと相関するという内生性の問題と、それに対処するBerry反転・操作変数法を学ぶ。今日の缶コーヒーの例も、次回さらに複雑にして再登場する。

参考文献

• Lancaster, Kelvin J. (1966) “A New Approach to Consumer Theory,” Journal of Political Economy, 74(2), 132-157. 論文
• 上武康亮・遠山祐太・若森直樹・渡辺安虎『実証ビジネス・エコノミクス』日本評論社、2025年、第2章。
• Train, Kenneth E. (2009) Discrete Choice Methods with Simulation, 2nd ed., Cambridge University Press, Chapters 2-3.
• McFadden, Daniel (1986) “The Choice Theory Approach to Market Research,” Marketing Science, 5(4), 275-297.
• McFadden, Daniel (2001) “Economic Choices,” American Economic Review, 91(3), 351-378.
• Guadagni, Peter M. and John D. C. Little (1983) “A Logit Model of Brand Choice Calibrated on Scanner Data,” Marketing Science, 2(3), 203-238.
• Manski, Charles F. (1977) “The Structure of Random Utility Models,” Theory and Decision, 8(3), 229-254. https://doi.org/10.1007/BF00133443