経済学分析でのモンテカルロシミュレーション：ダミーデータ生成編

この記事はCA Tech Lounge Advent Calendar 2024 13日目の記事です。昨日は@satousatouさんによる[chalice] dynamodbの変更をトリガーにAPI Gatewayのwebsocketで通知するでした。

こんにちは、@kod-yです。大学院で経済学、とりわけ実証産業組織論を専攻しています。本稿では、専攻分野でよく用いられるモンテカルロシミュレーションのうち、ダミーデータを生成するパートを実装例とともにご紹介します。

はじめに

まず、モンテカルロシミュレーションの概要と必要性についてご説明します。モンテカルロシミュレーションとは仮定したモデルに任意のパラメターを当ててシミュレーションし、ダミーデータを生成する手法です。他分野では目的変数の分布を確認する目的で行われるようですが、経済学では専ら実データを推定する前にモデルが正しくworkすることを確認する目的で行われます。

わざわざダミーデータまで作成して推定の練習をする原因は経済学モデルの複雑さにあります。経済学モデルは均衡状況の定式化のために需要側・供給側双方に多くの仮定を置きます。仮定から導出された式同士がどう関係しているかは可能な限り解析的に明らかにするものの、正しく定式化できているとは限りません。また、正しいモデルを実装したとしても、コードのバグで実は正しく推定できていないという可能性もあります。
そのため、実際に仮想的な経済主体をモデル通りに動かしてデータを生成し、正しく推定できることを確認します。これによって、実データが仮定した経済モデルに従う場合には推定が可能であると主張できます。

モデルの導入

He, Zheng, Belavina and Girotra (2020) の概要

本章では、本稿で議論するモデルを用意します。ここでは、He, Zheng, Belavina and Girotra (2020) のOnline Appendix 第7章にて説明されているシミュレーションを実装します。
He, Zheng, Belavina and Girotra (2020) の本文ではシェアサイクルプラットフォームにてユーザーがシェアサイクルを利用する要因としてステーションのネットワークが重要であることを主張するために、シェアサイクルの出発ステーションと到着ステーションの両方がわかっているデータ（route-levelデータ）を利用した離散選択モデルを推定しています。Online Appendixでは実モデルをより簡単化したモデルを用いて推定されたモデルが選択集合の内生性に対処できていることと、反実仮想分析がデータを地区ごとにaggregateしたデータの推定結果であっても妥当性を担保することを主張しています。

モデルの設定

先述の通り、モンテカルロシミュレーションでは何らかのtrue parameterをもとにモデルを動かしてダミーデータを生成します。He, Zheng, Belavina and Girotra (2020) ではモデルに従う個人が出発ステーションと到着ステーションの組み合わせ（ルート）を選択する状況を想定します。
具体的には、個人は以下のような効用関数を持つと仮定します。
$$
u_{iklt} = \beta_0 + \beta_1d_{kl} + \xi_{kl} + \epsilon_{iklt} \quad \mathrm{for}:kl \in C_t
$$
ここで $d_{kl}$ はルート $kl$ の距離、$\xi_{kl}$ はルート $kl$ の観察されない特性、$\epsilon_{iklt}$ はidiosyncratic shockを表しています。この効用関数を最大化するときに、$\epsilon_{iklt}$ が $i.i.d.$ に第一種極地分布に従うと仮定すると、個人 $i$ が時点 $t$ にルート $kl$ を選択する確率は以下のように表されます。
$$
p_{klt} = \frac{\exp(\beta_0 + \beta_1d_{kl} + \xi_{kl})}{1 + \sum_{k'l' \in C_t} \exp(\beta_0 + \beta_1 d_{k'l'} + \xi_{k'l'})}
$$

ダミーデータの生成

それでは、実際にダミーデータを生成します。
初めに、モデルに必要なパラメターのうち、期間やステーション数のような関心のないパラメターを設定します。

# number of customers
N <- 1000
# total time periods
T <- 1000
# number of stations
S <- 100

# map the stations on [0, 1] x [0, 1]
set.seed(2024)
data_station <- data.frame(
  station_id = 1:S,
  x = runif(S, 0, 1),
  y = runif(S, 0, 1)
)

次に、主要変数である $d_{kl}$ を設計します。Appendixではどのような距離を用いたのか言及がなかったので、stats::dist()で利用できるユークリッド距離を採用します。

# construct the distance matrix
dist_matrix <- as.matrix(dist(data_station[, c("x", "y")]))

更に、選択集合 $C_{it}$ の設計のために、bike availabilityとdock availabilityを設計します。Appendixでは簡単のために $[0.3, 0.7]$ の一様分布からドローしています。

# generate bike and dock availability
data_station <- data_station %>%
  mutate(
    availability_1_a = runif(S, 0.3, 0.7),
    availability_1_b = runif(S, 0.3, 0.7),
    availability_2_a = runif(S, 0.3, 0.7),
    availability_2_b = runif(S, 0.3, 0.7)
  ) %>%
  mutate(
    bike_availability_1 = pmin(availability_1_a, availability_1_b),
    dock_availability_1 = pmax(availability_1_a, availability_1_b),
    bike_availability_2 = pmin(availability_2_a, availability_2_b),
    dock_availability_2 = pmax(availability_2_a, availability_2_b)
  ) %>%
  select(-availability_1_a, -availability_1_b, -availability_2_a, -availability_2_b)

ここで2つずつドローしているのは、availabilityの一時的な揺れを許容するためです。新たに $\mathbb{P}{mix}$ を導入し、確率 $\mathbb{P}{mix}$ でbike_availability_1, dock_availability_1を、確率 $1 - \mathbb{P}_{mix}$ でbike_availability_2, dock_availability_2を採用します。

ここまででstation-levelのデータは作成できたので、これをroute-levelに拡張します。まず、station-levelデータをクロス結合し、ルート特有の観察されない特性を正規分布からドローします。

# cross join the station-level data
data_route <- cross_join(data_station, data_station)

# generate the unobserved characteristics
data_route <- data_route %>%
  mutate(
    unobserved = rnorm(S^2, 0, 1)
  )

結合後の変数を整理してから、ルート自体のavailability（$= C_t$）を定義します。

# define a route status based on p_mix and availability
data_BSS <- lapply(
  p_mix,
  function(p) {
    data_route_time %>%
      mutate(
        p_draw = runif(n(), 0, 1),
        bool_availablity = p_draw <= p
      ) %>%
      mutate(
        availability_status = case_when(
          bool_availablity & bike_availability_1.x < availability.x & dock_availability_1.y > availability.y ~ 1,
          !bool_availablity & bike_availability_2.x < availability.x & dock_availability_2.y > availability.y ~ 1,
          TRUE ~ 0
        )
      ) %>%
      select(
        station_id.x,
        station_id.y,
        time,
        distance,
        availability_status
      )
  }
)

次に、前章で定義した選択確率を求めます。

# derive a choice probability
# p_mix: c(0.5, 0.7, 0.9, 1)
# beta_1: c(-0.1, -0.25, -0.5, -0.75, -1, -2, -3, -5)
cp_true <- lapply(
  beta_1,
  function(beta) {
    lapply(
      data_BSS,
      function(data) {
        data %>%
          filter(
            availability_status == 1
          ) %>%
          mutate(
            utility = beta_0 + beta * distance,
            exp_utility = exp(utility)
          ) %>%
          group_by(time) %>%
          mutate(
            sum_exp_utility = 1 + sum(exp_utility)
          ) %>%
          ungroup() %>%
          mutate(
            choice_probability = exp_utility / sum_exp_utility
          ) %>%
          select(
            station_id.x,
            station_id.y,
            time,
            choice_probability
          )
      }
    )
  }
)

最後に、実際の選択をシミュレーションします。各 $t$ 期において選択確率に基づいたルート選択を $N$ 回行います。

# draw the choice
data_choice <- lapply(
  1:length(beta_1),
  function(i) {
    lapply(
      1:length(p_mix),
      function(j) {
        cp_true[[i]][[j]] %>%
          group_by(time) %>%
          sample_n(
            N, # size
            replace = TRUE, # replace is allowed
          )
      }
    )
  }
)

ここまでで、仮定したモデルに従ったときに生成されるダミーデータを得ることができました。

まとめ

本稿では、モンテカルロシミュレーションに必要なダミーデータを生成する手順をご紹介しました。それぞれの関数形の内容や推定方法は考えたいモデルによって異なりますが、大まかな手順は本モデルと同様になります。
また、紹介した実装につきましてはGitHubにてコードを公開する予定ですので必要に応じてご参照ください。

参考文献

• Pu He, Fanyin Zheng, Elena Belavina, Karan Girotra (2020) Customer Preference and Station Network in the London Bike-Share System. Management Science 67(3):1392-1412.