---
title: "課題8：staggered rolloutの効果検証と会員ランクRDD"
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(knitr)
library(fixest)
library(rdrobust)

set.seed(20268)

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)

# sunabの集計済み相対期間係数に対する、クラスタ頑健な同時Wald検定
sunab_pretrend_wald <- function(fit, periods, n_clusters) {
  beta <- fit$coefficients
  beta_names <- names(beta)
  V_raw <- vcov(fit)[beta_names, beta_names, drop = FALSE]
  X_raw <- model.matrix(fit)[, beta_names, drop = FALSE]
  period_name <- sub("::.*$", "", rownames(aggregate(fit, "period"))[1])
  targets <- paste0(period_name, "::", periods)
  targets <- targets[vapply(
    targets,
    function(z) any(startsWith(beta_names, paste0(z, ":cohort::"))),
    logical(1)
  )]

  A <- matrix(0, nrow = length(targets), ncol = length(beta),
              dimnames = list(targets, beta_names))
  for (j in seq_along(targets)) {
    idx <- startsWith(beta_names, paste0(targets[j], ":cohort::"))
    shares <- colSums(sign(X_raw[, idx, drop = FALSE]))
    A[j, idx] <- shares / sum(shares)
  }

  b <- drop(A %*% beta)
  V <- A %*% V_raw %*% t(A)
  eig <- eigen(V, symmetric = TRUE)
  keep <- eig$values > max(eig$values) * 1e-8
  V_inv <- eig$vectors[, keep, drop = FALSE] %*%
    diag(1 / eig$values[keep], nrow = sum(keep)) %*%
    t(eig$vectors[, keep, drop = FALSE])
  df1 <- sum(keep)
  F_stat <- drop(t(b) %*% V_inv %*% b) / df1

  data.frame(
    tested_periods = paste(periods, collapse = ", "),
    F_statistic = F_stat,
    df1 = df1,
    df2 = n_clusters - 1,
    p_value = pf(F_stat, df1, n_clusters - 1, lower.tail = FALSE)
  )
}
```

## この課題について

第8回の講義では、staggered adoption（施策の導入時期が単位ごとにバラバラ）のもとで、素朴な静的TWFE（two-way fixed effects）推定がどれほど深刻に歪みうるか、そして`fixest::sunab`でSun-Abraham型のevent studyを推定する方法をシミュレーションで確認した。また、会員ランクのような閾値を利用したRDDの基本も学んだ。

この課題では、講義ノートと**同じ構造・異なる設定**のDGPを自分の手で構築し、静的TWFEとSun-Abraham型推定を比較した上で、会員ランクRDDの推定を独立に行う。

- **Part 1**：小売チェーンの**店舗**（都道府県ではなく個別店舗）における新レジシステムの時差導入DGPを自作する（配点40点）。
- **Part 2**：(i) 前後比較・単純比較、(ii) 静的TWFE、(iii) sunabによるevent study、(iv) pre-trend確認、(v) 会員ランクRDD、(vi) 経営含意3行、を行う（配点60点）。

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

::: {.callout-important}
## この課題のゴール

- staggered adoption・動的異質効果というDGPを自分の手で書けるようになる。
- 静的TWFEが真の効果の符号やスケールをどう見誤るかを、自分のシミュレーションで確認する。
- `fixest::sunab`を使ったevent studyの実装と、その出力の読み方（係数名の対応、`ref.p`の意味）を身につける。
- `rdrobust`を使った基本的なRDD推定と、走行変数のヒストグラムをbunchingの診断材料として読む作法を身につける。
:::

## Part 1：店舗レベルstaggered rolloutのシミュレーション設計（配点40点）

### 舞台設定

ある全国スーパーマーケットチェーンが、**新しいセルフレジシステム**を店舗ごとに時差を持って導入したとする。店舗を $i = 1, \dots, G$（$G = 47$、講義ノートと同じ店舗数を採用するが、ここでは「都道府県」ではなく「店舗」という単位に読み替える）、月を $t = 1, \dots, T$（$T = 36$、3年分の月次データ）とする。

新レジシステムの効果は、次のような構造を持つとする。

- 導入時期は店舗ごとにバラバラ（staggered）。導入時点を $t = 6$ から $t = 30$ の範囲でランダムに割り当てる。
- 約20%の店舗は、観測期間中ずっと未導入（never-treated）のままとする。
- 効果は動的：導入後の経過月数 $e = t - \text{adopt}_i$ に対して線形に成長する。ただし成長の傾き $\theta_i$ は店舗ごとに異なり、**早期導入店舗ほど効果の伸びが大きく**（$\theta_{\text{early}} = 0.18$）、**後期導入店舗ほど効果の伸びが小さい**（$\theta_{\text{late}} = -0.10$、講義ノートの値から少し変えてある）。傾きは導入時点に応じて線形補間で決まる：

$$
\theta_i = \theta_{\text{early}} + \frac{\text{adopt}_i - \text{adopt}_{\text{lo}}}{\text{adopt}_{\text{hi}} - \text{adopt}_{\text{lo}}} \left(\theta_{\text{late}} - \theta_{\text{early}}\right)
$$

- 店舗固定効果 $\mu_i \sim N(0, 1)$、月固定効果 $\lambda_t \sim N(0, 0.4^2)$、共通トレンド $0.012 \times t$、誤差項は標準偏差0.25の正規分布。

アウトカム $y_{it}$（例えば1人あたりレジ通過時間の短縮指数、あるいは顧客満足度指数など、具体的な解釈は自由に設定してよい）は次のモデルに従う。

$$
y_{it} = \mu_i + \lambda_t + 0.012 t + \theta_i \cdot e_{it} \cdot D_{it} + \varepsilon_{it}
$$

ここで $D_{it} = 1\{t \geq \text{adopt}_i\}$、$e_{it} = t - \text{adopt}_i$（$D_{it}=1$のときのみ意味を持つ）。

::: {.callout-tip}
## ヒント：DGPを書く順番

講義ノートのRシミュ2と同じ手順でよい。

1. 店舗数`G_a`・観測月数`T_a`・never-treated比率を決める。
2. `sample()`で導入時点を店舗ごとにランダムに割り当てる（never-treated店舗には観測期間を大きく超える値、例えば`T_a + 1000`を割り当てるとよい）。
3. 導入時点から線形補間で店舗固有の傾き`theta_i`を計算する。
4. 店舗固定効果・月固定効果を生成する。
5. `expand.grid()`で店舗×月のパネルを作り、`treated`・`e`（経過月数）・`true_effect`・`y`を計算する。
:::

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

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

- 店舗数 $G = 47$、観測月数 $T = 36$
- never-treated比率：20%
- 導入時点の範囲：6〜30ヶ月目
- $\theta_{\text{early}} = 0.18$、$\theta_{\text{late}} = -0.10$
- 店舗固定効果の標準偏差：1.0、月固定効果の標準偏差：0.4
- 共通トレンド係数：0.012、誤差項の標準偏差：0.25

```{r}
#| eval: false
G_a <- 47
T_a <- 36
never_frac_a <- 0.20
adopt_lo_a <- 6
adopt_hi_a <- 30

n_never_a <- round(G_a * never_frac_a)
n_treated_a <- G_a - n_never_a

adopt_time_treated_a <- sample(adopt_lo_a:adopt_hi_a, size = n_treated_a, replace = TRUE)
adopt_time_a <- c(adopt_time_treated_a, rep(T_a + 1000, n_never_a))
adopt_time_a <- ____________________  # 店舗にランダムに割り当て直す（sample()を使う）

theta_early_a <- 0.18
theta_late_a  <- -0.10

theta_i_a <- ifelse(
  adopt_time_a > T_a,
  0,
  ____________________  # 線形補間の式（本文の数式を参考に）
)

store_fe_a <- rnorm(G_a, mean = 0, sd = 1.0)
month_fe_a <- rnorm(T_a, mean = 0, sd = 0.4)

rollout_df <- expand.grid(store = 1:G_a, month = 1:T_a) %>%
  mutate(
    adopt_t = adopt_time_a[store],
    theta_store = theta_i_a[store],
    treated = as.numeric(month >= adopt_t),
    e = ifelse(treated == 1, month - adopt_t, NA_real_),
    true_effect = ifelse(treated == 1, theta_store * e, 0),
    y = store_fe_a[store] + month_fe_a[month] + 0.012 * month +
      true_effect + rnorm(n(), mean = 0, sd = 0.25)
  )

head(rollout_df)
```

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

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

(a) 早期導入店舗（導入時点が最も早い店舗）1店と、後期導入店舗（導入時点が最も遅い店舗）1店の、月次アウトカムの推移を折れ線グラフで比較する図。
(b) 全店舗の導入時点のヒストグラム（横軸：導入月、never-treated店舗は除く）。
(c) 処置セル（`treated == 1`）における`true_effect`を、経過月数`e`ごとに平均して折れ線グラフで描く（真の動的効果の形状を可視化する）。

::: {.callout-tip}
## ヒント：導入時点の早い店舗・遅い店舗の抽出

```r
early_store <- rollout_df %>% filter(adopt_t <= T_a) %>%
  filter(adopt_t == min(adopt_t)) %>% pull(store) %>% unique() %>% first()
late_store <- rollout_df %>% filter(adopt_t <= T_a) %>%
  filter(adopt_t == max(adopt_t)) %>% pull(store) %>% unique() %>% first()
```

のように、`adopt_t`の最小値・最大値を持つ店舗を1つずつ選べばよい。
:::

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

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

- never-treated店舗の数と、処置群店舗の導入時点の平均・最小・最大。
- 処置セル（`treated == 1`）における`true_effect`の平均・標準偏差（これが後で「真のATT」として使われる）。
- 早期コホート（導入時点が中央値より早い店舗群）と後期コホート（中央値より遅い店舗群）それぞれの`theta_store`の平均値（コホート間の異質性を数値で確認する）。

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

### 設問2-1（10点）：前後比較・単純比較

処置群店舗について、導入前（`treated == 0`の期間）と導入後（`treated == 1`の期間）の`y`の平均を単純に比較せよ。また、**最終月**における処置群店舗と対照群（never-treated）店舗の平均`y`をそのまま比較せよ。

これら2つの「素朴な」比較の値を、Part 1で確認した真のATT（`true_effect`の平均）と並べて表にまとめ、どちらの比較がより真値に近いか、あるいはどちらも大きくずれているかをコメントせよ。

::: {.callout-tip}
## ヒント：なぜこれらの比較は不十分か

前後比較は、共通トレンド（`0.012 * month`）や月固定効果を処置効果に混ぜる。最終月の群間比較は店舗の水準差を混ぜるうえ、最終月の効果と「全処置セルの平均ATT」という異なる対象を比べている。数値が偶然真値に近いかどうかではなく、各比較が何を混ぜ、どの対象を測っているかを説明するのが狙いである。次の静的TWFEも万能な正解ではなく、動学的・コホート別異質性のため別の歪みを持つ。
:::

### 設問2-2（12点）：静的TWFE推定

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

```r
static_fit_a <- feols(y ~ treated | store + month,
                      data = rollout_df, vcov = ~store)
```

推定された係数を、Part 1で確認した真のATT（`true_effect`の平均）と比較し、**符号が一致しているか、大きさはどの程度ズレているか**を明示的にコメントせよ。

::: {.callout-tip}
## ヒント：符号のズレに気づいたら

もし静的TWFEの符号が真のATTと逆になっていたら、それは実装ミスではなく、講義で学んだ「禁じられた比較」（すでに処置され、かつ動学的に効果が成長し続けている店舗を対照に使ってしまう）が起きている可能性が高い。慌てず、講義ノートの説明を振り返りながら、なぜこの現象が起きるのかを自分の言葉で書けるようにしておこう。
:::

### 設問2-3（12点）：sunabによるevent studyの推定

`fixest::sunab()`を使って、コホート×相対期間の異質性に頑健なSun-Abraham型event studyを推定せよ。

```r
sunab_fit_a <- feols(y ~ sunab(adopt_t, month) | store + month,
                     data = rollout_df, vcov = ~store)
```

以下を行え。

(a) `iplot(sunab_fit_a)`でevent studyの図を描く。
(b) 静的TWFEの推定値と、sunabの全体ATT集計値（`summary(sunab_fit_a, agg = "att")`で得られる）を並べた比較表を作り、真のATTとあわせて3つを比較する図（棒グラフ）を描く。

::: {.callout-tip}
## ヒント：sunabの係数名

`coef(sunab_fit_a)`で係数名を確認すると、`month::-3`や`month::5`のような形式で「処置からの相対期間」に対応する係数が並んでいるはずである。負の数字が処置前（pre-period）、正の数字（および0）が処置後（post-period）に対応する。
:::

### 設問2-4（8点）：pre-trendの記述と同時検定

sunabの推定結果について、次の2つを行え。

1. 直近4期の処置前係数（$e=-5,-4,-3,-2$）を取り出し、平均絶対値を計算する。
2. セットアップチャンクで用意した`sunab_pretrend_wald()`を使い、この4係数が同時にゼロかを、店舗クラスタ頑健なWald検定で確かめる。

図、平均絶対値、Wald検定のp値をまとめて解釈せよ。「棄却できないことは並行トレンドの証明ではない」ことと、pre-trend検定は小さな違反に対して検出力が低いことも明記すること。

::: {.callout-tip}
## ヒント：処置前係数の抽出

```r
#| eval: false
pre_periods_a <- -5:-2
period_table_a <- aggregate(sunab_fit_a, "period")
pre_names_a <- paste0("month::", pre_periods_a)
pre_coefs_a <- period_table_a[pre_names_a, "Estimate"]
mean(abs(pre_coefs_a), na.rm = TRUE)
sunab_pretrend_wald(sunab_fit_a, pre_periods_a, n_clusters = G_a)
```
:::

### 設問2-5（別課題・13点）：会員ランクRDD

ここからは別のデータセットを使う独立した設問である。ある小売チェーンの会員制度で、**年間購入額が5万円以上になるとゴールド会員になり、翌年の購買額が変化する**という状況を考える。次のDGPでデータを生成せよ（このDGPには走行変数の操作＝bunchingは含めない）。

```{r}
#| eval: false
n_rdd_a <- 3000
cutoff_a <- 5.0

x_rdd_a <- rgamma(n_rdd_a, shape = 3.0, scale = 1.8)
x_rdd_a <- pmin(pmax(x_rdd_a, 0.1), 20)

true_jump_a <- 1.2

baseline_a <- function(x) {
  3.0 + 0.6 * (x - cutoff_a) - 0.03 * (x - cutoff_a)^2
}

gold_a <- as.numeric(x_rdd_a >= cutoff_a)
y_rdd_a <- baseline_a(x_rdd_a) + true_jump_a * gold_a + rnorm(n_rdd_a, mean = 0, sd = 0.8)

rdd_df_a <- data.frame(
  annual_purchase = x_rdd_a,
  gold_status = gold_a,
  next_year_purchase = y_rdd_a
)
```

以下を行え。

(a) `rdrobust::rdrobust(y, x, c = cutoff)`でジャンプを推定し、`Conventional`の点推定値と`Robust`の95%信頼区間を報告して、真値1.2万円と比較せよ。
(b) `rdrobust::rdplot(y, x, c = cutoff)`でビン分割散布図を描け。
(c) 走行変数`annual_purchase`のヒストグラムを、カットオフ付近に注目して描き、カットオフ直下・直上のカウント比を計算せよ。この標本で目立つbunchingが見えないかを診断し、単純な比が1に近くても「操作なし」やRDDの妥当性の証明にはならないことを述べよ。
(d) この推定値を因果効果として読むために必要な仮定を二つ挙げ、推定対象がカットオフ近傍の局所効果であることを説明せよ。

::: {.callout-tip}
## ヒント：カットオフ直下・直上のカウント

```r
#| eval: false
count_below <- sum(rdd_df_a$annual_purchase >= cutoff_a - 0.5 & rdd_df_a$annual_purchase < cutoff_a)
count_above <- sum(rdd_df_a$annual_purchase >= cutoff_a & rdd_df_a$annual_purchase < cutoff_a + 0.5)
count_above / count_below
```

操作なしのDGPなので、この比はおおむね1に近い値になるはずである。もし3倍・4倍のように大きくずれていたら、DGPの実装を見直すこと。ただし、これはビン幅に左右される粗い診断であり、McCrary型の密度検定そのものではない。
:::

### 設問2-6（経営含意・5点）：経営・政策への含意を3行で書け

あなたは、この小売チェーンのデジタル施策担当マネージャーである。上司から「新レジシステムの効果検証、静的な前後比較でざっくり出しておいて」と言われた。

この依頼に対して、**3行程度**で返信メモを書け。静的な前後比較（あるいは静的TWFE）のどこに注意が必要で、代わりにどのような分析を提案するか、簡潔に述べること。

::: {.callout-tip}
## ヒント：メモに何を書くべきか

(1) staggered rolloutの場合、静的な比較は符号すら誤ることがあるという事実を一言で説明する、(2) 実務でどの分析（sunabのようなevent study）を使うべきか明言する、(3) 次のロールアウトでは導入順序をランダム化する提案を添える、という3点を意識するとよい。
:::

## 提出物

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