--- title: "課題1：小売パネルデータで価格弾力性を推定する" subtitle: "計量経済学II" author: "Kei Ikegami" lang: ja format: html: toc: true toc-depth: 3 toc-location: left number-sections: false embed-resources: true code-copy: true theme: cosmo fig-width: 8 fig-height: 4.8 execute: warning: false message: false --- ```{r} #| include: false library(ggplot2) library(dplyr) library(tidyr) library(fixest) set.seed(20261) theme_lecture <- function() { fam <- if (Sys.info()[["sysname"]] == "Darwin") "Hiragino Sans" else "" theme_minimal(base_size = 13, base_family = fam) } theme_set(theme_lecture()) knitr::opts_chunk$set(fig.align = "center", dpi = 144) options(scipen = 8) ``` ## この課題について 第1回の講義では、「値上げすると販売数量が増える？」というシミュレーションを通じて、価格の内生性（繁忙期には価格も需要も同時に上がる）がOLSによる弾力性推定をどう歪めるかを見た。この課題では、その内容を**小売チェーンの店舗パネルデータ**（複数店舗×複数週）という、実務でより典型的な設定に拡張し、計量経済学Iで学んだ**操作変数法（2SLS）**を使ってこの内生性に対処する練習をする。 やることは大きく2つである。 - **Part 1**：自分でデータ生成過程（DGP）を書き、20店舗×52週分の小売パネルデータをシミュレーションで作る。 - **Part 2**：そのデータを使って、(i) OLSで弾力性を推定してバイアスを確認し、(ii) 操作変数法（2SLS）で真の弾力性を回収し、(iii) モンテカルロ・シミュレーションで両者の標本分布を比較する。 配点は合計100点。手を動かせば2〜4時間程度で終わる分量である。各設問にヒントcalloutを用意しているので、詰まったら遠慮なく開いてほしい。 ::: {.callout-important} ## この課題のゴール - パネルデータにおける価格の内生性（需要ショックが価格と数量の両方に影響する）を、自分の手で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} = \mu + \eta \ln p_{it} + \gamma_t + \mu_i + \varepsilon_{it} $$ - $\eta = -1.5$：真の価格弾力性（既知の値としてシミュレーションに使う。推定する対象）。 - $\gamma_t$：季節的な需要ショック（週ごとに共通、繁忙期・閑散期の周期的な変動）。 - $\mu_i$：店舗固定効果（店舗ごとの需要水準の違い、立地や客層の違いなど）。 - $\varepsilon_{it}$：観測されない週次のノイズ。 価格は次のように決まるとする。 $$ \ln p_{it} = \pi_0 + \pi_1 c_{it} + \pi_2 \gamma_t + \nu_{it} $$ - $c_{it}$：**仕入コストのシフター**（原材料費指数など、店舗・週ごとに変動する費用側の変数）。需要とは無関係に決まる、価格に影響を与える変数。 - $\gamma_t$：Part 1のモデルで説明した季節的な需要ショックと**同じ変数**。つまり店舗（の価格設定担当者）は、繁忙期には価格を強気に設定する。 - $\nu_{it}$：価格設定における観測されないノイズ。 このモデルの構造を図にすると、次のような関係になる。$\gamma_t$（需要ショック）が $\ln p_{it}$ と $\ln q_{it}$ の**両方**に影響を与えている点に注目してほしい。これが価格の内生性の源泉であり、$c_{it}$（仕入コスト）は $\ln p_{it}$ には影響するが $\ln q_{it}$ には（価格を経由する以外は）影響しないという点が、操作変数として使える理由である。 ::: {.callout-tip} ## ヒント：DGPを書く順番 以下の順番でコードを組み立てると迷いにくい。 1. 店舗ID・週IDのベクトルを作る（`rep()`, `rep_len()`などを使う。`n_stores * n_weeks`行のデータフレームを想定）。 2. 店舗固定効果 $\mu_i$ を20店舗分生成する（`rnorm`）。 3. 季節ショック $\gamma_t$ を52週分生成する。単純な正弦波（`sin(2*pi*week/52)`）に週次ノイズを乗せるとよい。 4. 仕入コストシフター $c_{it}$ を生成する（`rnorm`、店舗×週で変動）。 5. 価格 $\ln p_{it}$ を、$c_{it}$ と $\gamma_t$ の線形結合として生成する。 6. 需要 $\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 ```{r} #| eval: false 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つの図を作成せよ。 (a) 価格と需要（対数）の散布図（回帰直線付き）。全店舗・全週をプールしたもの。 (b) ある1店舗（`store_id == 1`など）について、週ごとの価格の推移と需要の推移を、それぞれ折れ線グラフで示す（2つの図でも、`patchwork`で並べても良い）。 (c) 価格のヒストグラム、需要のヒストグラム。 ::: {.callout-tip} ## ヒント：折れ線グラフ `dplyr::filter(panel_df, store_id == 1)` で1店舗分に絞り込んでから、`ggplot(aes(x = week_id, y = log_price)) + geom_line()` のように描けばよい。価格と需要が同じような山谷のタイミングで動いていることが確認できるはずである（これが繁忙期に価格も需要も一緒に動くというDGPの反映である）。 ::: ### 設問1-3（15点）：記述統計 `dplyr::summarise()` や `knitr::kable()` を使って、次を表にまとめよ。 - 価格・需要それぞれの平均、標準偏差、最小値、最大値（全店舗・全週プール）。 - 価格と需要（levelでも対数でもよい）の単純な相関係数。 - 価格と仕入コストシフターの相関係数、価格と需要ショックの相関係数（この2つの数値が、後で2SLSがなぜ必要かを理解する伏線になる）。 ## 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$ と比較せよ。バイアスの大きさと方向（上方か下方か）についてコメントすること。 ::: {.callout-tip} ## ヒント：なぜバイアスが出るか 価格は需要ショック $\gamma_t$ に正に反応して決まる（$\pi_2 = 0.25 > 0$）。需要ショックは需要そのものにも正の影響を与える。つまり「価格が高い週＝需要ショックが大きい週」であり、需要ショックの効果が価格の係数に混ざり込む。真の弾力性は負なので、この混入は係数を**真の値よりも大きい方向（0に近い、あるいは正）**にバイアスさせる。 ::: ### 設問2-2（20点）：操作変数法（2SLS）による弾力性の回収 仕入コストシフター $c_{it}$ は、価格には影響するが（$\pi_1 = 0.15 \neq 0$）、需要には価格を経由する以外の経路で影響しない、という設定になっている。つまり $c_{it}$ は価格 $\ln p_{it}$ の**操作変数（instrumental variable, IV）**として使える。 計量経済学Iで学んだ2SLSの考え方を復習する。 1. **第1段階**：$\ln p_{it}$ を $c_{it}$（と固定効果）に回帰し、予測値 $\widehat{\ln p_{it}}$ を得る。 2. **第2段階**：$\ln q_{it}$ を $\widehat{\ln p_{it}}$（と固定効果）に回帰する。 `fixest::feols()` のIV構文を使うと、この2段階を自動でやってくれる。構文は次の形になる。 ```{r} #| eval: false # IV構文： y ~ 外生説明変数 | 固定効果 | 内生変数 ~ 操作変数 iv_fit <- feols(log_demand ~ 1 | store_id | log_price ~ cost_shifter, data = panel_df) summary(iv_fit) ``` 以下を行え。 (a) 上記の2SLS推定を実行し、$\hat\eta_{IV}$ を得て、OLS推定値・真の値と並べた比較表を作る。 (b) 第1段階の回帰（`log_price ~ cost_shifter`、固定効果込み）を単独で実行し、操作変数の関連性（instrument relevance）を、係数の大きさと有意性で確認せよ。 (c) （発展・任意）第1段階のF統計量を確認し、「弱操作変数」の心配がないかコメントせよ。 ::: {.callout-tip} ## ヒント：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推定値の平均は真の値にどれだけ近いか。 - 両者の標準偏差（ばらつき）を比較せよ。一般に操作変数法は、内生性を除去する代わりに推定量のばらつきが大きくなることが知られているが、操作変数（仕入コストシフター）の説明力が強い場合はこの分散拡大は小さくなる。手元の結果がどちらの方向に出たか、その理由（第1段階の関連性の強さ）とあわせてコメントせよ。 ::: {.callout-tip} ## ヒント：モンテカルロのコード構造 計量経済学Iのモンテカルロ演習と同じ構造でよい。 ```{r} #| eval: false 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くらいと出たので、この数字を使って値上げの影響を試算してほしい」と言われた。 この依頼に対して、**3行程度**で返信メモを書け。何が問題で、どの数字を代わりに使うべきで、なぜそう言えるのかを簡潔に述べること。 ::: {.callout-tip} ## ヒント：メモに何を書くべきか 単に「OLSはダメでIVが正しい」ではなく、(1) なぜOLSの-0.6が過小評価（絶対値として小さい）なのか一言で説明する、(2) 実務でどの数字を使うべきか明言する、(3) 可能であれば「他にもコストシフターのような操作変数を使えるデータがないか」という前向きな提案を添える、という3点を意識すると実務的なメモになる。 ::: ## 提出物 - 本ファイル（`assignment1.qmd`）をレンダリングしたHTMLファイル。 - 全てのコードチャンクが上から順に実行可能であること（`eval: false` のチャンクは自分で `eval: true` に直すか、コードをコピーして実行すること）。