課題1：小売パネルデータで価格弾力性を推定する

計量経済学II

作者

Kei Ikegami

この課題について

第1回の講義では、広告の例で予測と介入後のシミュレーションの違いを学び、価格実験から推定した需要を使って値付けを提案した。この課題では、価格をランダムに割り当てた実験とは異なり、店舗が需要の動きを見て値付けをしている小売チェーンの店舗パネルデータ（複数店舗×複数週）を扱う。計量経済学Iで学んだ操作変数法（2SLS）を使って価格の内生性に対処し、推定した弾力性を値上げの影響試算につなげる練習である。

やることは大きく2つである。

Part 1：自分でデータ生成過程（DGP）を書き、20店舗×52週分の小売パネルデータをシミュレーションで作る。
Part 2：そのデータを使って、(i) OLSと操作変数法（2SLS）で弾力性を推定し、(ii) モンテカルロ・シミュレーションで両者の標本分布を比較し、(iii) 推定値から値上げ時の数量変化を試算する。

配点は合計100点。手を動かせば2〜4時間程度で終わる分量である。各設問にヒントcalloutを用意しているので、詰まったら遠慮なく開いてほしい。

この課題のゴール

パネルデータにおける価格の内生性（需要ショックが価格と数量の両方に影響する）を、自分の手でDGPとして書けるようになる。
log-log OLSによる弾力性推定がどの方向にバイアスを持つかを、シミュレーションで確認する。
コストシフター（費用側の変数）を操作変数として使った2SLS推定を、fixest::feolsのIV構文で実装できるようになる。
モンテカルロ・シミュレーションによって、OLSとIVの推定量の標本分布の違い（バイアスの有無）を可視化できるようになる。

Part 1：小売パネルデータのシミュレーション設計（配点40点）

モデルの設定

20店舗、各店舗52週（1年分）のデータを考える。店舗を \(i = 1, \dots, 20\)、週を \(t = 1, \dots, 52\) とする。

対数需要（販売数量）は次のモデルに従うとする。

\[ \ln q_{it} = \ln 150 + \eta \ln(p_{it}/500) + \gamma_t + \mu_i + \varepsilon_{it} \]

\(\eta = -1.5\)：真の価格弾力性（既知の値としてシミュレーションに使う。推定する対象）。
\(\gamma_t\)：季節的な需要ショック（週ごとに共通、繁忙期・閑散期の周期的な変動）。
\(\mu_i\)：店舗固定効果（店舗ごとの需要水準の違い、立地や客層の違いなど）。
\(\varepsilon_{it}\)：観測されない週次のノイズ。

価格が500円、\(\gamma_t=\mu_i=0\) のとき、数量の中央値が週150個となるように切片を設定している。

価格は次のように決まるとする。

\[ \ln p_{it} = \pi_0 + \pi_1 c_{it} + \pi_2 \gamma_t + \nu_{it} \]

\(c_{it}\)：仕入コストのシフター（原材料費指数など、店舗・週ごとに変動する費用側の変数）。需要とは無関係に決まる、価格に影響を与える変数。
\(\gamma_t\)：Part 1のモデルで説明した季節的な需要ショックと同じ変数。つまり店舗（の価格設定担当者）は、繁忙期には価格を強気に設定する。
\(\nu_{it}\)：価格設定における観測されないノイズ。

\(\pi_0=\ln 500\) とし、需要ショックは \(\gamma_t=0.3\sin(2\pi t/52)+w_t\) とする。乱数部分は

\[ \begin{aligned} \mu_i&\sim N(0,0.3^2), & w_t&\sim N(0,0.25^2), & c_{it}&\sim N(0,1),\\ \varepsilon_{it}&\sim N(0,0.1^2), & \nu_{it}&\sim N(0,0.05^2) \end{aligned} \]

に従い、これらの乱数は互いに独立に生成する。同じ店舗の行は \(\mu_i\) を、同じ週の行は \(\gamma_t\) を共有する。

\(\gamma_t\)（需要ショック）が \(\ln p_{it}\) と \(\ln q_{it}\) の両方に影響を与えている点に注目してほしい。これを回帰から落とすと価格に内生性が生じる。一方、\(c_{it}\)（仕入コスト）は需要側の要因とは独立で、数量には価格を経由する以外に影響しない。この設定と価格への関連性が、操作変数として使える理由である。

真の弾力性や demand_shock を知っているのは、ここでは自分でデータを作るからである。真値はデータ生成と推定結果の検証に使い、Part 2の推定では需要ショックを観測できないと考える。最後の値上げの試算には、真値ではなくデータから得た推定値を使う。

ヒント：DGPを書く順番

以下の順番でコードを組み立てると迷いにくい。

店舗ID・週IDのベクトルを作る（rep(), rep_len()などを使う。n_stores * n_weeks行のデータフレームを想定）。
店舗固定効果 \(\mu_i\) を20店舗分生成する（rnorm）。
季節ショック \(\gamma_t\) を52週分生成する。単純な正弦波（sin(2*pi*week/52)）に週次ノイズを乗せるとよい。
仕入コストシフター \(c_{it}\) を生成する（rnorm、店舗×週で変動）。
価格 \(\ln p_{it}\) を、\(c_{it}\) と \(\gamma_t\) の線形結合として生成する。
需要 \(\ln q_{it}\) を、真の弾力性 \(\eta=-1.5\) を使って \(\ln p_{it}\)、\(\gamma_t\)、\(\mu_i\) の線形結合として生成する。

\(\gamma_t\) は店舗によらず「週」だけで決まる変数なので、店舗IDでループを回すのではなく、週ID変数から作った \(\gamma_t\) を、各店舗の行に対応させる（rep()で店舗数だけ繰り返すか、週IDをキーにして結合する）とよい。

設問1-1（10点）：パラメータとDGPコードの実装

以下の骨格を埋めて、20店舗×52週のパネルデータを生成せよ。パラメータは次の値を使うこと。

真の価格弾力性 \(\eta = -1.5\)
価格が仕入コストに反応する強さ \(\pi_1 = 0.15\)
価格が需要ショックに反応する強さ \(\pi_2 = 0.25\)（この値が0でない、というのが内生性の源泉）
店舗固定効果の標準偏差：0.3
需要側ノイズ \(\varepsilon_{it}\) の標準偏差：0.1
価格側ノイズ \(\nu_{it}\) の標準偏差：0.05

n_stores <- 20
n_weeks <- 52
N <- n_stores * n_weeks

store_id <- rep(1:n_stores, each = n_weeks)
week_id  <- rep(1:n_weeks, times = n_stores)

# 店舗固定効果（20店舗、sd=0.3）
store_fe <- ____________________
store_fe_long <- store_fe[store_id]

# 季節的な需要ショック（52週分。正弦波 + ノイズ）
season_wave <- 0.3 * sin(2 * pi * (1:n_weeks) / 52)
demand_shock_weekly <- season_wave + rnorm(n_weeks, mean = 0, sd = 0.25)
demand_shock <- demand_shock_weekly[week_id]

# 仕入コストシフター（店舗×週で変動する外生変数）
cost_shifter <- rnorm(N, mean = 0, sd = 1)

# 真のパラメータ
true_elasticity <- -1.5
pi1_cost <- 0.15
pi2_demand <- 0.25

# 価格の生成
log_price <- log(500) + pi1_cost * cost_shifter + pi2_demand * demand_shock +
  rnorm(N, mean = 0, sd = 0.05)
price <- exp(log_price)

# 需要の生成
log_demand <- ____________________
demand <- exp(log_demand)

panel_df <- data.frame(
  store_id = store_id,
  week_id = week_id,
  cost_shifter = cost_shifter,
  demand_shock = demand_shock,
  price = price,
  log_price = log_price,
  demand = demand,
  log_demand = log_demand
)

head(panel_df)

設問1-2（15点）：可視化

以下の3つの図を作成せよ。

価格と需要（対数）の散布図（回帰直線付き）。全店舗・全週をプールしたもの。
ある1店舗（store_id == 1など）について、週ごとの価格の推移と需要の推移を、それぞれ折れ線グラフで示す（2つの図でも、patchworkで並べても良い）。
価格のヒストグラム、需要のヒストグラム。

ヒント：折れ線グラフ

dplyr::filter(panel_df, store_id == 1) で1店舗分に絞り込んでから、ggplot(aes(x = week_id, y = log_price)) + geom_line() のように描けばよい。需要ショックは価格と数量の両方を押し上げるが、仕入コストなども毎週変わるので、実際の価格と数量の山谷が常に一致するわけではない。

設問1-3（15点）：記述統計

dplyr::summarise() や knitr::kable() を使って、次を表にまとめよ。

価格・需要それぞれの平均、標準偏差、最小値、最大値（全店舗・全週プール）。
価格と需要（levelでも対数でもよい）の単純な相関係数。
価格と仕入コストシフターの相関係数、価格と需要ショックの相関係数（価格を動かす2つの要因を確認する。相関だけで操作変数の妥当性を証明できるわけではない）。

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

設問2-1（15点）：log-log OLSによる弾力性推定

Part 1で作ったデータを使い、次の回帰を実行せよ。

\[ \ln q_{it} = a + \eta \ln p_{it} + u_{it} \]

店舗固定効果 \(\mu_i\) を明示的にコントロールする場合としない場合の両方を推定し（lm() または fixest::feols()、固定効果ありは feols(log_demand ~ log_price | store_id, data = panel_df) の形式）、推定された \(\hat\eta\) を真の値 \(\eta = -1.5\) と比較せよ。真値からのずれの大きさと方向を述べ、DGPから予想されるバイアスの方向と対応させること。1回の推定のずれと、繰り返し標本を生成したときのバイアスは区別する。

ヒント：なぜバイアスが出るか

価格は需要ショック \(\gamma_t\) に正に反応して決まる（\(\pi_2 = 0.25 > 0\)）。需要ショックは需要そのものにも正の影響を与える。価格が高い週ほど需要ショックも大きい傾向があるため、これを回帰から落とすと、その効果が価格の係数に混ざり込む。真の弾力性は負なので、この混入は係数を真の値よりも大きい方向（0に近い、あるいは正）にバイアスさせる。

なお、このDGPでは需要ショックが全店舗で共通なので、店舗固定効果に加えて週固定効果も入れれば \(\gamma_t\) を取り除ける（| store_id + week_id）。IVだけが対処法ではない。以下では操作変数法の練習として、店舗固定効果のみを入れたOLSと2SLSを比較する。

設問2-2（20点）：操作変数法（2SLS）による弾力性の回収

仕入コストシフター \(c_{it}\) は、価格には影響するが（\(\pi_1 = 0.15 \neq 0\)）、需要には価格を経由する以外の経路で影響しない、という設定になっている。つまり \(c_{it}\) は価格 \(\ln p_{it}\) の操作変数（instrumental variable, IV）として使える。

計量経済学Iで学んだ2SLSの考え方を復習する。

第1段階：\(\ln p_{it}\) を \(c_{it}\)（と固定効果）に回帰し、予測値 \(\widehat{\ln p_{it}}\) を得る。
第2段階：\(\ln q_{it}\) を \(\widehat{\ln p_{it}}\)（と固定効果）に回帰する。

fixest::feols() のIV構文を使うと、この2段階を自動でやってくれる。構文は次の形になる。

# IV構文： y ~ 外生説明変数 | 固定効果 | 内生変数 ~ 操作変数
iv_fit <- feols(log_demand ~ 1 | store_id | log_price ~ cost_shifter,
                data = panel_df)
summary(iv_fit, vcov = ~week_id)

以下を行え。

上記の2SLS推定を実行し、\(\hat\eta_{IV}\) を得て、OLS推定値・真の値と並べた比較表を作る。
第1段階の回帰（log_price ~ cost_shifter、店舗固定効果込み）を単独で実行し、操作変数の関連性（instrument relevance）を、係数の大きさと有意性で確認せよ。さらに、関連性以外に必要な条件がこのDGPでどう満たされているか説明せよ。
（発展・任意、配点外）第1段階のF統計量を確認し、「弱操作変数」の心配がないかコメントせよ。

同じ週の店舗は需要ショックを共有するので、(a)(b)の標準誤差は summary(推定結果, vcov = ~week_id) として週単位でクラスター化する。

ヒント：feolsのIV構文でよく間違えるところ

feols() のIV構文は y ~ 外生変数 | 固定効果 | 内生変数 ~ 操作変数 という3段のパイプ区切りである。外生の説明変数が無い場合は y ~ 1 | 固定効果 | 内生変数 ~ 操作変数 のように 1 を置く。固定効果が無いモデルなら真ん中のパイプを省略して y ~ 1 | 内生変数 ~ 操作変数 と書く。summary(iv_fit) はデフォルトでは第2段階の結果だけを表示するので、第1段階を見るには summary(iv_fit, stage = 1) を使うか、設問(b)のように第1段階の回帰を別途 feols() で単独に実行するとよい。

設問2-3（15点）：モンテカルロ・シミュレーションによる標本分布の比較

Part 1のDGPを200回繰り返してデータを生成し直し、そのたびにOLS推定値とIV推定値を計算せよ。両方とも店舗固定効果のみを入れる。得られた200個のOLS推定値、200個のIV推定値をヒストグラム（または密度プロット）で重ねて描き、真の値 \(\eta = -1.5\) に縦線を引いて比較せよ。

以下を確認・報告すること。

OLS推定値の平均は真の値からどれだけ乖離しているか（バイアス）。
IV推定値の平均は真の値にどれだけ近いか。
両者の標準偏差（ばらつき）とRMSE（二乗誤差の平均の平方根）を比較せよ。偏ったOLSとIVの比較では、標準偏差の大小だけで推定の良さを判断できない。真値からのずれとばらつきを合わせてコメントせよ。

ヒント：モンテカルロのコード構造

計量経済学Iのモンテカルロ演習と同じ構造でよい。

n_mc <- 200
ols_estimates <- numeric(n_mc)
iv_estimates  <- numeric(n_mc)

for (m in seq_len(n_mc)) {
  # 1. Part1と同じDGPで、ただしseedを変えてデータを再生成する
  # 2. OLS推定 -> ols_estimates[m] に格納
  # 3. IV推定  -> iv_estimates[m] に格納
}

ループの中で set.seed() を呼ぶ場合は、set.seed(1000 + m) のように毎回異なる値にすること（同じ値を使うと200回とも同じデータになってしまう）。1チャンクの実行時間の目安は30秒以内、ノート全体では数分以内に収まるよう、n_mc <- 200 程度に留めておくこと。

設問2-4（10点）：値付け担当者への3行メモ

あなたはこの小売チェーンの価格戦略チームに所属するアナリストである。上司（値付け担当のシニアマネージャー）から「うちのデータでざっくり価格弾力性を見て、OLSで-0.6くらいと出たので、この数字を使って値上げの影響を試算してほしい」と言われた。

同じ店舗・同じ週の需要条件を保って価格を10%上げる場合を考え、設問2-2で推定したIV弾力性を使って数量の変化率を試算せよ。その結果を含めて、3行程度で返信メモを書け。OLSの数字の問題、試算した数量変化、試算を信頼するための条件を簡潔に述べること。

ヒント：メモに何を書くべきか

一定弾力性の需要式では、価格を1.1倍にしたときの数量変化率は \(100\{1.1^{\hat\eta_{IV}}-1\}\)% である。データ生成に使った \(-1.5\) を代入するのではなく、自分の推定値を使う。この数量比では、需要条件と需要側ノイズの分布が変わらなければ、対数から水準に戻す際の補正項も分子・分母で相殺される。

メモでは、(1) 需要ショックの混入によってOLSが値上げによる数量減を小さく見積もりうること、(2) IV推定値から計算した数量変化、(3) 仕入コストが需要に直接影響せず需要側の要因と独立であることに加え、値上げ後も需要式が変わらないという前提を意識するとよい。

提出物

本ファイル（assignment1.qmd）をレンダリングしたHTMLファイル。
全てのコードチャンクが上から順に実行可能であること（eval: false のチャンクは自分で eval: true に直すか、コードをコピーして実行すること）。