1 はじめに

2022年6・7月号の「実証ビジネス・エコノミクス」の第8回「価格戦略をダイナミックに考える：シングルエージェント動学モデルの推定［応用編］」に付随するRコードの紹介になります。

1.1 留意事項

本連載の内容、およびサンプルコード等の資料は、情報の提供のみを目的としていますので、運用につきましては十分にご確認をいただき、お客様ご自身の責任とご判断によって行ってください。これらの情報の運用結果により損害等が生じた場合でも、日本評論社および著者はいかなる責任を負うことはできませんので、ご留意ください。

1.2 謝辞

今回のプログラムの作成に際して、島本幸典さん(東京大学大学院経済学研究科)にご尽力頂きました。この場を借りて御礼申し上げます。

1.3 全体の流れ

下準備：パッケージの導入、データの読み込み
記述統計の確認
二段階推定方法の実践
反実仮想分析

2 下準備

2.1 Rに関する下準備

# ワークスペースを初期化
rm(list = ls())

# パッケージを読み込む
require(tidyverse)
require(skimr)
require(fixest)
require(numDeriv)
require(kableExtra)

2.2 データの読み込み

第8回では、第7回で生成した自動車の買い替え行動に関するデータを引き続き用いる。

# 第7回で生成したデータを読み込む。
data_gen <- 
  readr::read_csv('Chap8_data.csv')

## Rows: 120000 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (5): consumer, period, action, price, mileage
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

data_gen %>% head(3)

## # A tibble: 3 × 5
##   consumer period action price mileage
##      <dbl>  <dbl>  <dbl> <dbl>   <dbl>
## 1        1    481      0  2400      95
## 2        1    482      1  2100      95
## 3        1    483      0  2300       0

前回同様、状態変数とその番号の関係を示すデータフレームを作成する。

## Stateの作成

# 価格の状態変数
price_states <- seq(2000, 2500, by = 100)

# 走行距離の状態変数
mileage_states <- seq(0, 100, by = 5)

# 価格の状態変数の数
num_price_states <- length(price_states)

# 走行距離の状態変数の数
num_mileage_states <- length(mileage_states)

# 状態変数は価格と走行距離の状態変数のペア
# 従って状態変数の数は価格の状態変数の数と走行距離の状態変数の数の積となる
num_states <- num_price_states * num_mileage_states

# 価格、走行距離の状態変数の組み合わせ(p,m)を1つのデータフレームで表す
state_df <- 
  dplyr::tibble(
    # stateにidを番号付ける
    state_id = 1:num_states,
    # 順番は (p,m) = (2000,0),(2100,0),...,(2500,0),(2000,5),(2100,5),...
    price_id = rep(1:num_price_states, times = num_mileage_states),
    mileage_id = rep(1:num_mileage_states, each = num_price_states),
    price = rep(price_states, times = num_mileage_states),
    mileage = rep(mileage_states, each = num_price_states)
  )

# 下3行を表示
state_df %>% tail(3)

## # A tibble: 3 × 5
##   state_id price_id mileage_id price mileage
##      <int>    <int>      <int> <dbl>   <dbl>
## 1      124        4         21  2300     100
## 2      125        5         21  2400     100
## 3      126        6         21  2500     100

# 分析のためにデータを加工する
data_gen <- 
  data_gen %>% 
  # 状態変数を表した列を追加
  dplyr::left_join(state_df, by = c('price', 'mileage')) %>% 
  dplyr::group_by(consumer) %>% 
  # 遷移行列の推定で使うため、ラグ変数を追加
  dplyr::mutate(lag_price_id = lag(price_id),
                lag_mileage_id = lag(mileage_id),
                lag_action = lag(action)) %>% 
  dplyr::ungroup()

3 記述統計の確認

生成したデータの記述統計を確認する。

# 生成したデータの要約統計
data_gen %>% 
  dplyr::select(price, mileage, action) %>%
  skimr::skim() %>% 
  skimr::yank("numeric") %>% 
  dplyr::select(skim_variable, mean, sd, p0, p100)

Variable type: numeric

skim_variable	mean	sd	p0	p100
price	2323.71	152.03	2000	2500
mileage	33.24	23.83	0	100
action	0.03	0.17	0	1

4 二段階推定方法の実践

# 推定に必要なパラメタの設定

# 時間割引率
beta <- 0.99

# オイラー定数
Euler_const <- - digamma(1)

# 消費者は車を購入するかどうか決めるため選択肢の数は2つ
num_choice <- 2

4.1 Step 1: CCP の推定、及び State transition の推定

4.1.1 State transition の推定

第7回と同様に遷移行列をデータから推定する。遷移行列を推定するコードは前回と全く同じである。まず、走行距離の遷移行列を推定する。

# それぞれの確率が実現した観察の数を数える
num_cond_obs_mileage <- 
  data_gen %>% 
  # 1期目は推定に使えないため落とす
  dplyr::filter(period != min(data_gen$period)) %>% 
  # t期の走行距離、t+1期の走行距離、t期の購買ごとにグループ化して、観察数を数える
  dplyr::group_by(lag_mileage_id, mileage_id, lag_action) %>% 
  dplyr::summarise(num_cond_obs = n(),
                   .groups = 'drop') %>% 
  # 確率ごとに名前を割り当てる
  dplyr::mutate(
    cond_obs_mileage = case_when(
      # 1 - kappa_1 - kappa_2 の場合
      (
        (lag_action == 0 &
           between(lag_mileage_id, 1, 20) &
           (lag_mileage_id == mileage_id)) |
          (lag_action == 1 & 
             mileage_id == 1)
      ) ~ 'cond_obs_mileage1',
      # kappa_1 の場合
      (
        (lag_action == 0 &
           between(lag_mileage_id, 1, 19) &
           (lag_mileage_id == mileage_id - 1)) |
          (lag_action == 1 & 
             mileage_id == 2)
      ) ~ 'cond_obs_mileage2',
      # kappa_2 の場合
      (
        (lag_action == 0 &
           between(lag_mileage_id, 1, 19) &
           (lag_mileage_id == mileage_id - 2)) |
          (lag_action == 1 & 
             mileage_id == 3)
      ) ~ 'cond_obs_mileage3',
      # kappa_1 + kappa_2 の場合
      (
        lag_action == 0 &
          lag_mileage_id == 20 &
          mileage_id == 21
      ) ~ 'cond_obs_mileage4',
      TRUE ~ 'other'
    )) %>% 
  # 'other' は推定には使わないため落とす
  dplyr::filter(cond_obs_mileage != 'other') %>% 
  # 確率ごとにグループ化し、再度、観察の数を数える
  dplyr::group_by(cond_obs_mileage) %>% 
  dplyr::summarise(num_cond_obs = as.numeric(sum(num_cond_obs)),
                   .groups = 'drop') %>% 
  dplyr::select(num_cond_obs) %>% 
  as.matrix()

# 最尤法の解析解により推定値を求める
kappa_est <- c()
kappa_est[1] <- 
  (num_cond_obs_mileage[2] * 
     (num_cond_obs_mileage[2] + num_cond_obs_mileage[3] + num_cond_obs_mileage[4])) /
  ((num_cond_obs_mileage[2] + num_cond_obs_mileage[3]) * 
     (num_cond_obs_mileage[1] + num_cond_obs_mileage[2] + 
        num_cond_obs_mileage[3] + num_cond_obs_mileage[4]))
kappa_est[2] <- 
  (num_cond_obs_mileage[3] * 
     (num_cond_obs_mileage[2] + num_cond_obs_mileage[3] + num_cond_obs_mileage[4])) /
  ((num_cond_obs_mileage[2] + num_cond_obs_mileage[3]) * 
     (num_cond_obs_mileage[1] + num_cond_obs_mileage[2] + 
        num_cond_obs_mileage[3] + num_cond_obs_mileage[4]))

次に、価格の遷移行列を推定する。

# それぞれの確率が実現した観察の数を数える
num_cond_obs_price <- 
  data_gen %>% 
  # 1期目は推定に使えないため落とす
  dplyr::filter(period != min(data_gen$period)) %>% 
  # t期の価格、t+1期の価格ごとにグループ化して、観察数を数える
  dplyr::group_by(lag_price_id, price_id) %>% 
  dplyr::summarise(num_cond_obs = n(),
                   .groups = 'drop') %>% 
  # 観察数を行列（num_price_states行の正方行列）に変換
  # price_id (t+1期の価格) を横に広げる
  tidyr::pivot_wider(names_from = "price_id",
                     values_from = "num_cond_obs") %>%
  dplyr::select(!lag_price_id) %>% 
  as.matrix()

# 最尤法の解析解により推定値を求める
lambda_est_mat <- 
  num_cond_obs_price / rowSums(num_cond_obs_price)
lambda_est_mat

##               1          2          3          4         5          6
## [1,] 0.10273973 0.10388128 0.19931507 0.19874429 0.1921233 0.20319635
## [2,] 0.10074627 0.09916327 0.19222071 0.20952058 0.1991180 0.19923112
## [3,] 0.10191150 0.10086410 0.29714585 0.19780047 0.2038230 0.09845509
## [4,] 0.09703103 0.10031847 0.19359975 0.30331827 0.2055681 0.10016437
## [5,] 0.04924386 0.04994318 0.09886652 0.09840030 0.5007722 0.20277397
## [6,] 0.05057165 0.05123799 0.09826752 0.09931963 0.2030932 0.49751000

lambda_est <- as.vector(t(lambda_est_mat))[c(-1,-8,-15,-22,-29,-36)]

走行距離、価格の遷移行列を生成する関数を作成し、推定された遷移行列を取得する。

# パラメタを所与として走行距離の遷移行列を出力する関数を作成
gen_mileage_trans <- function(kappa){
  # 走行距離が1,2段階上がる確率をパラメタとする
  kappa_1 <- kappa[1]
  kappa_2 <- kappa[2]
  # 購買しなかった場合の遷移行列を作成
  mileage_trans_mat_hat_not_buy <-
    matrix(0, ncol = num_mileage_states, nrow = num_mileage_states)
  for (i in 1:num_mileage_states) {
    for (j in 1:num_mileage_states) {
      if (i == j){
        mileage_trans_mat_hat_not_buy[i,j] <- 1 - kappa_1 - kappa_2
      } else if (i == j - 1) {
        mileage_trans_mat_hat_not_buy[i,j] <- kappa_1
      } else if (i == j - 2){
        mileage_trans_mat_hat_not_buy[i,j] <- kappa_2
      }
    }
  }
  mileage_trans_mat_hat_not_buy[num_mileage_states - 1, num_mileage_states] <- 
    kappa_1 + kappa_2
  mileage_trans_mat_hat_not_buy[num_mileage_states, num_mileage_states] <- 1
  # 購買した場合の遷移行列を作成
  # 購入した期では m=0 となるため次の期のmileageはそこから決まることに注意
  mileage_trans_mat_hat_buy <-
    matrix(1, nrow = num_mileage_states, ncol = 1) %*%
    mileage_trans_mat_hat_not_buy[1,]
  # 3次元のarrayとして出力
  return(array(c(mileage_trans_mat_hat_not_buy, 
                 mileage_trans_mat_hat_buy), 
               dim=c(num_mileage_states,num_mileage_states,num_choice)))
}

# パラメタを所与として価格の遷移行列を出力する関数を作成
gen_price_trans <- function(lambda){
  lambda_11 <- 1 - lambda[1] - lambda[2] - lambda[3] - lambda[4] - lambda[5]
  lambda_22 <- 1 - lambda[6] - lambda[7] - lambda[8] - lambda[9] - lambda[10]
  lambda_33 <- 1 - lambda[11] - lambda[12] - lambda[13] - lambda[14] - lambda[15]
  lambda_44 <- 1 - lambda[16] - lambda[17] - lambda[18] - lambda[19] - lambda[20]
  lambda_55 <- 1 - lambda[21] - lambda[22] - lambda[23] - lambda[24] - lambda[25]
  lambda_66 <- 1 - lambda[26] - lambda[27] - lambda[28] - lambda[29] - lambda[30]
  price_trans_mat_hat <- 
    c(lambda_11, lambda[1], lambda[2], lambda[3], lambda[4], lambda[5],
      lambda[6], lambda_22, lambda[7], lambda[8], lambda[9], lambda[10],
      lambda[11], lambda[12], lambda_33, lambda[13], lambda[14], lambda[15],
      lambda[16], lambda[17], lambda[18], lambda_44, lambda[19], lambda[20],
      lambda[21], lambda[22], lambda[23], lambda[24], lambda_55, lambda[25],
      lambda[26], lambda[27], lambda[28], lambda[29], lambda[30], lambda_66) %>% 
    matrix(ncol = num_price_states, nrow = num_price_states, byrow=T)
  return(price_trans_mat_hat)
}

# 推定された遷移行列を取得（本誌のG(i)に対応）
# 配列の第三次元が購買を意味しており、i=0,1という順番で並んでいる。
G <- array(c(gen_mileage_trans(kappa_est)[,,1] %x% gen_price_trans(lambda_est), 
             gen_mileage_trans(kappa_est)[,,2] %x% gen_price_trans(lambda_est)),
           dim=c(num_states, num_states, num_choice))

4.1.2 CCPの推定

今回はロジットモデルによる近似で、CCPを推定する。

# 単純なロジットモデルを推定する
logit_model <- fixest::feglm(action ~ price + price^2 + mileage + mileage^2,
                             data_gen, family = binomial("logit"))
# logit_model <- feglm(action ~ price + price^2 + 
#                        mileage + mileage^2 +
#                        price*mileage + (price*mileage)^2
#                      , data_gen, family = binomial("logit"))
# logit_model_glm <-  glm(action ~ price + price^2 + mileage + mileage^2, 
#                         data_gen, family = binomial("logit"))

# 推定値に基づいてCCPを予測する。列はi=0,1という順番。
CCP_1st <- 
  cbind(
    1 - predict(logit_model, state_df),
    predict(logit_model, state_df)
  )

4.2 Step 2: パラメタの推定

以下のコードでは、二段階目の方法として紹介のあった (1) 行列形式によるインバージョン、と (2) 有限依存性 (finite dependence) アプローチの２つの手法を実装する。

4.2.1 行列形式によるインバージョン

# 前回同様、状態変数、コントロール変数毎の今期の効用を返す関数を作成
flow_utility <- function(theta, state_df){
  theta_c <- theta[1]
  theta_p <- theta[2]
  U <- 
    cbind(
      # その期における車を購入しない場合の効用
      U_not_buy = - theta_c * state_df$mileage, 
      
      # その期における車を購入する場合の効用
      U_buy = - theta_p * state_df$price
    ) 
  return(U)
}

# 前回同様、行列の(i,j)要素を出力する関数を作成
mat_ij <- Vectorize(
  function(i,j,mat) {mat[i,j]},
  vectorize.args = c("i", "j"))

# 行列インバージョンで期待価値関数を導出し、CCPを出力する関数を作成
policy_operator_mat_inv <- function(theta, CCP, beta, G, state_df){
  # 今期の効用を計算（本誌のu(i)に対応）
  U <- flow_utility(theta, state_df)
  # psiを計算（本誌のpsi(i)に対応）
  num_states <- nrow(state_df)
  psi <- Euler_const * matrix(1, num_states, num_choice) - log(CCP)
  # 行列計算により期待価値関数ベクトルを計算
  V <- 
    solve(diag(num_states) - 
            beta * (CCP[,1] %*% matrix(1, 1, num_states) * G[,,1] + 
                      CCP[,2] %*% matrix(1, 1, num_states) * G[,,2])) %*% 
    rowSums(CCP * (U + psi))
  # 選択肢ごとの価値関数を計算
  CV <- U + beta * cbind(G[,,1] %*% V, G[,,2] %*% V)
  # CCPを更新
  CCP <- exp(CV) / rowSums(exp(CV))
  return(CCP)
}

likelihood_fun <- function(theta, CCP, df, beta, G, state_df, policy_operator){
  # CCP operaterにより、CCPを更新する
  # 今回は行列インバージョンと有限依存性(finite dependence)アプローチのどちらか
  CCP <- policy_operator(theta, CCP, beta, G, state_df)
  # 疑似尤度を計算
  obj <- sum(log(mat_ij(df$state_id, df$action + 1, CCP)))
  return(obj)
}

# パラメタの真の値
theta_true <- c(theta_c = 0.004, theta_p = 0.003)

start_time_mat_inv <- proc.time()

# 最適化
mat_inv_opt_mat_inv <- optim(theta_true, likelihood_fun,
                             CCP = CCP_1st, df = data_gen, 
                             beta = beta, G = G, state_df = state_df, 
                             policy_operator = policy_operator_mat_inv,
                             control = list(fnscale = -1), 
                             method = "Nelder-Mead")

end_time_mat_inv <- proc.time()
run_time_mat_inv <- (end_time_mat_inv - start_time_mat_inv)[[3]]
print(stringr::str_c("Runtime: ", round(run_time_mat_inv,2)))

## [1] "Runtime: 15.08"

# 結果の確認
theta_mat_inv <- mat_inv_opt_mat_inv$par
theta_mat_inv

##     theta_c     theta_p 
## 0.003884738 0.002971597

# 前回と同様、標準誤差を推定する
hessian <- numDeriv::hessian(func = likelihood_fun, x = theta_mat_inv,
                             CCP = CCP_1st, df = data_gen, 
                             beta = beta, G = G, state_df = state_df, 
                             policy_operator = policy_operator_mat_inv)
theta_se_mat_inv <- sqrt(diag(solve(-hessian)))
dplyr::tibble(theta_mat_inv, theta_se_mat_inv)

## # A tibble: 2 × 2
##   theta_mat_inv theta_se_mat_inv
##           <dbl>            <dbl>
## 1       0.00388        0.0000846
## 2       0.00297        0.0000358

4.2.2 有限依存性 (finite dependence) アプローチ

# 行列インバージョンで期待価値関数を導出し、CCPを出力する関数を作成
policy_operator_finite_dep <- function(theta, CCP, beta, G, state_df){
  # 今期の効用を計算（本誌のu(i)に対応）
  U <- flow_utility(theta, state_df)
  # Finite dependenceに基づき、選択ごとの価値関数の差を計算
  CV_dif <- 
    U[,2] - U[,1] + 
    beta * (G[,,2] %*% (-log(CCP[,2])) - G[,,1] %*% (-log(CCP[,2])))
  # CCPを更新
  prob_buy <- exp(CV_dif) / (1 + exp(CV_dif))
  CCP <- cbind(1-prob_buy, prob_buy)
  return(CCP)
}

start_time_finite_dep <- proc.time()

# 最適化
finite_dep_opt <- optim(theta_true, likelihood_fun,
                        CCP = CCP_1st, df = data_gen, 
                        beta = beta, G = G, state_df = state_df, 
                        policy_operator = policy_operator_finite_dep,
                        control = list(fnscale = -1), 
                        method = "Nelder-Mead")

end_time_finite_dep <- proc.time()
run_time_finite_dep <- (end_time_finite_dep - start_time_finite_dep)[[3]]
print(stringr::str_c("Runtime: ", round(run_time_finite_dep,2)))

## [1] "Runtime: 9.79"

# 結果の確認
theta_finite_dep <- finite_dep_opt$par
theta_finite_dep

##     theta_c     theta_p 
## 0.003821963 0.002903782

# 標準誤差を推定する
hessian <- numDeriv::hessian(func = likelihood_fun, x = theta_finite_dep,
                             CCP = CCP_1st, df = data_gen, 
                             beta = beta, G = G, state_df = state_df, 
                             policy_operator = policy_operator_finite_dep)
theta_se_finite_dep <- sqrt(diag(solve(-hessian)))
dplyr::tibble(theta_finite_dep, theta_se_finite_dep)

## # A tibble: 2 × 2
##   theta_finite_dep theta_se_finite_dep
##              <dbl>               <dbl>
## 1          0.00382           0.000793 
## 2          0.00290           0.0000219

5 推定方法の比較

5.1 第7回の復習：入れ子不動点アルゴリズム

以下は、本誌第8回第3節のノーテーションに合わせたコードになっている点に注意されたい。そのため、前回（第7回）と比べると若干コードに変更点があるが、本質的には全く同じである。当然、推定値も同じ結果が得られている。

contraction <- 
  function(theta, beta, G, state_df) {
    # 価値関数の初期値
    num_states <- nrow(state_df)
    V_old <- matrix(0, nrow = num_states, ncol = 1)
    # パラメタより今期の効用を計算
    U <- flow_utility(theta, state_df)
    
    # 価値関数の差の初期値
    diff <- 1000
    # 縮小写像の誤差範囲
    tol_level <- 1.0e-10
    
    while (diff > tol_level) {
      # 価値関数を計算
      V_new <- log(rowSums(exp(U + beta * cbind(G[,,1] %*% V_old, G[,,2] %*% V_old)))) + Euler_const
      # 価値関数の更新による差を評価
      diff <- sum(abs(V_new-V_old))
      # 価値関数を次のループに渡す(価値関数の更新)
      V_old <- V_new
    }
    # 得られた事前の価値関数を出力
    return(V_old)
  }

policy_operator_nfxp <- 
  function(theta, beta, G, state_df) {
    # パラメタより今期の効用を計算
    U <- flow_utility(theta, state_df)
    # 縮小写像アルゴリズムで得られた事前の価値関数をVとする
    V <- contraction(theta, beta, G, state_df)
    # 選択肢ごとの価値関数を計算
    CV <- U + beta * cbind(G[,,1] %*% V, G[,,2] %*% V)
    # CCPを計算
    CCP <- exp(CV) / rowSums(exp(CV))
    return(CCP)
  }

likelihood_fun_nfxp <- function(theta, df, beta, G, state_df){
  # CCP operaterにより、
  CCP <- policy_operator_nfxp(theta, beta, G, state_df)
  # 疑似尤度を計算
  obj <- sum(log(mat_ij(df$state_id, df$action + 1, CCP)))
  return(obj)
}

start_time_nfxp <- proc.time()
# 最適化
nfxp_opt <- optim(theta_true, likelihood_fun_nfxp,
                  df = data_gen, 
                  beta = beta, G = G, state_df = state_df, 
                  control = list(fnscale = -1), 
                  method = "Nelder-Mead")

end_time_nfxp <- proc.time()
run_time_nfxp <- (end_time_nfxp - start_time_nfxp)[[3]]
print(stringr::str_c("Runtime: ", round(run_time_nfxp,2)))

## [1] "Runtime: 37.08"

theta_nfxp <- nfxp_opt$par
theta_nfxp

##     theta_c     theta_p 
## 0.003887129 0.002971684

# 標準誤差を推定する
hessian <- numDeriv::hessian(func = likelihood_fun_nfxp, x = theta_finite_dep,
                             df = data_gen, 
                             beta = beta, G = G, state_df = state_df)
theta_se_nfxp <- sqrt(diag(solve(-hessian)))
dplyr::tibble(theta_nfxp, theta_se_nfxp)

## # A tibble: 2 × 2
##   theta_nfxp theta_se_nfxp
##        <dbl>         <dbl>
## 1    0.00389     0.0000854
## 2    0.00297     0.0000348

5.2 推定方法の比較表

dplyr::tibble(
  algorithm = c('NFXP', 'Matrix', 'Finite', 'True'),
  theta_c = c(theta_nfxp[1], theta_mat_inv[1], theta_finite_dep[1], theta_true[1]),
  theta_se_c = c(theta_se_nfxp[1], theta_se_mat_inv[1], theta_se_finite_dep[1], NA),
  theta_p = c(theta_nfxp[2], theta_mat_inv[2], theta_finite_dep[2], theta_true[2]),
  theta_se_p = c(theta_se_nfxp[2], theta_se_mat_inv[2], theta_se_finite_dep[2], NA),
  run_time = c(run_time_nfxp, run_time_mat_inv, run_time_finite_dep, NA)
) %>%
  kableExtra::kable(format = "html",
                    booktabs = TRUE,
                    caption = "推定方法の比較",
                    digits = c(0,6,6,6,6,2), 
                    format.args = list(scientific = FALSE)) %>%
  kableExtra::kable_styling()

推定方法の比較
algorithm	theta_c	theta_se_c	theta_p	theta_se_p	run_time
NFXP	0.003887	0.000085	0.002972	0.000035	37.08
Matrix	0.003885	0.000085	0.002972	0.000036	15.08
Finite	0.003822	0.000793	0.002904	0.000022	9.79
True	0.004000	NA	0.003000	NA	NA

6 反実仮想分析1: Every Day Low Price

6.1 EDLPの効果

以下の二つのシナリオについてそれぞれCCPを計算し、結果を比較する。

推定した際と全く同じ状況（以降、ベースラインと呼ぶ）
価格が最も低い価格に永続的に固定される状況

# 以降のシミュレーションで得られたCCPを要素とするリストを作成
CCP_list <- list()

# シナリオ1) 推定した際と全く同じ状況
CCP_list[['Baseline']] <- 
  policy_operator_nfxp(theta_nfxp, beta, G, state_df)

# シナリオ1における、各価格の購買確率を計算
result_df_edlp <- 
  data_gen %>% 
  # データよりそれぞれの状態変数が観察された数を数える
  dplyr::group_by(state_id, price_id, price) %>% 
  dplyr::summarise(num_obs = n(),
                   .groups = 'drop') %>% 
  dplyr::arrange(state_id) %>% 
  # シナリオ1で得られたCCPを新たな列に加える
  dplyr::mutate(prob_buy_baseline = CCP_list[['Baseline']][,2]) %>% 
  # 走行距離について、購買確率の加重平均を取る（走行距離の頻度について積分する）
  dplyr::group_by(price_id, price) %>% 
  dplyr::summarize(prob_buy_baseline = weighted.mean(prob_buy_baseline, num_obs), 
                   .groups = 'drop')

# シナリオ2) 価格が最も低い価格に永続的に固定される状況

# このシナリオでは価格の変動がないため、走行距離についてのみの遷移行列がこのモデルの遷移行列となる
G_fixed_price <- gen_mileage_trans(kappa_est)

# 価格を固定した場合のCCPを計算
# 後ほどの需要と収入の計算で使うため、価格を固定した場合のCCPをすべての価格に対して計算
for (fixed_price in price_states) {
  # 価格が固定された場合の、状態変数を示すデータフレームを作成
  state_df_fixed_price <- 
    state_df %>% 
    dplyr::filter(price == fixed_price) %>% 
    dplyr::arrange(mileage_id)
  
  # 遷移行列と状態変数を示すデータフレームが変わっていることに注意してCCPを計算
  CCP_list[[stringr::str_c('edlp', fixed_price)]] <- 
    policy_operator_nfxp(theta = theta_nfxp, beta = beta, 
                         G = G_fixed_price, state_df = state_df_fixed_price)
}

# 200万円に価格を固定した場合の購買確率を計算
prob_buy_edlp2000 <- 
  data_gen %>% 
  # データよりそれぞれの状態変数が観察された数を数える
  dplyr::group_by(mileage_id, mileage) %>% 
  dplyr::summarise(num_obs = n(),
                   .groups = 'drop') %>% 
  dplyr::arrange(mileage_id) %>% 
  # 価格が固定された場合のCCPを新たな列に加える
  dplyr::mutate(prob_buy = CCP_list[['edlp2000']][,2]) %>% 
  # 走行距離について、購買確率の加重平均を取る（走行距離の頻度について積分する）
  dplyr::summarize(prob_buy = weighted.mean(prob_buy, num_obs)) %>% 
  as.matrix()

# 価格ごとの購買確率を図示する
result_df_edlp %>% 
  ggplot2::ggplot(aes(x = price, y = prob_buy_baseline)) + 
  ggplot2::geom_bar(stat = "identity") +
  ggplot2::geom_hline(yintercept = prob_buy_edlp2000,
                      color = "blue", size = 1.0) +
  ggplot2::annotate("text", 
                    x = 2500, y = prob_buy_edlp2000 + 0.002, 
                    color= "blue", 
                    label = str_c("EDLP: ", round(prob_buy_edlp2000,5))) +
  ggplot2::theme_minimal()

7 反実仮想分析2: 永続的・一時的な値下げ

以下の３つのシナリオを比較することで、値下げの効果を分析する。

ベースライン
10万円の永続的な値下げ
10万円の一時的な値下げ

7.1 値下げが起きた時のCCPの計算

まず、それぞれのシナリオにおけるCCPを計算する。

# シナリオ2) 10万円の永続的な値下げ
# priceを100下げた、state_dfを用意し、CCPを再計算する
state_df_discount <- 
  state_df %>% 
  dplyr::mutate(price = price - 100)

CCP_list[['Permanent']] <- 
  policy_operator_nfxp(theta_nfxp, beta, G, state_df_discount)

# シナリオ3) 10万円の一時的な値下げ

# 値下げが起きたときの今期の効用を計算
U_discount <- flow_utility(theta_nfxp, state_df_discount)
# 元の価格設定における事前の価値関数を計算
V <- contraction(theta_nfxp, beta, G, state_df)
# 選択肢ごとの価値関数を計算
CV_temporary <- U_discount + beta * cbind(G[,,1] %*% V, G[,,2] %*% V)
# CCPを計算
CCP_list[['Temporary']] <-exp(CV_temporary) / rowSums(exp(CV_temporary))

# 三つのシナリオについて、元々の価格が220万円の場合の、走行距離ごとの購入確率を図示
state_df %>% 
  # 得られたCCPを一つのdataframeにまとめる
  dplyr::mutate(ProbBaseline = CCP_list[['Baseline']][,2],
                ProbPermanent = CCP_list[['Permanent']][,2],
                ProbTemporary = CCP_list[['Temporary']][,2]) %>% 
  # 220万円での購買確率に絞る
  dplyr::filter(price == 2200) %>% 
  dplyr::select(mileage, ProbBaseline, 
                ProbPermanent, ProbTemporary) %>% 
  tidyr::pivot_longer(cols = - mileage) %>% 
  dplyr::rename(scenario = name, prob_buy = value) %>% 
  ggplot2::ggplot(aes(x = mileage, y = prob_buy, color = scenario)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() +
  ggplot2::theme_minimal()

7.2 消費者の分布のシミュレーション（値下げ）

最初の期における、状態変数ごとの消費者の分布を所与とした時に、条件付き選択確率（CCP）と遷移行列を使えば、時間を通じた消費者の分布の遷移をシミュレートできる。また、その分布とCCPを使うことで、各期の購買割合や収入がわかる。これらの計算を用いて、値下げがデータ上の購買割合に与える影響や、需要と収入への影響を観察する。

# データから状態変数毎の消費者の分布を計算
# state_idの順番に並んだ、(1 * num_states) の行ベクトル
consumer_dist_obs <- 
  data_gen %>% 
  dplyr::group_by(state_id, price_id, mileage_id) %>% 
  dplyr::summarise(num_obs = n(),
                   .groups = 'drop') %>% 
  dplyr::mutate(consumer_dist_obs = num_obs/sum(num_obs)) %>% 
  dplyr::arrange(state_id) %>% 
  dplyr::select(consumer_dist_obs) %>% 
  as.matrix() %>% 
  t()

# 購買を考慮した遷移行列を作成
discount_scenerio_vec <- c('Baseline','Permanent','Temporary')

# CCPを考慮した遷移行列を要素とするリストを作成
G_CCP_list <- list()

# それぞれの値下げシナリオに応じた、購買を考慮した遷移行列を作成
for (scenario in discount_scenerio_vec) {
  G_CCP_list[[scenario]] <- 
    (CCP_list[[scenario]][,1] %*% matrix(1, 1, num_states)) * G[,,1] +
    (CCP_list[[scenario]][,2] %*% matrix(1, 1, num_states)) * G[,,2]
}

# 需要と収入の計算に使うため、EDLPを行ったときの購買を考慮した遷移行列も作成
for (scenario in stringr::str_c('edlp', price_states)) {
  G_CCP_list[[scenario]] <- 
    (CCP_list[[scenario]][,1] %*% matrix(1, 1, nrow(G_fixed_price[,,1]))) * G_fixed_price[,,1] +
    (CCP_list[[scenario]][,2] %*% matrix(1, 1, nrow(G_fixed_price[,,1]))) * G_fixed_price[,,2]
}

# シミュレーションする消費者数と期間を定義
num_consumer_sim <- 1000
num_period_sim <- 20

# 値下げによるシミュレーション結果を表すデータフレームを前もって作成
discount_sim_df <- 
  dplyr::tibble(period = 1:num_period_sim)

# それぞれのシナリオについてシミュレーションを行う
for (scenario in discount_scenerio_vec) {
  
  # 消費者の分布の推移を表す行列を前もって作成
  consumer_dist_sim <- matrix(0, num_period_sim, num_states)
  # 初期の分布をデータ上の分布とする
  consumer_dist_sim[1,] <- consumer_dist_obs
  
  # 購買割合、需要、収入の空の列を作成
  discount_sim_df <-
    discount_sim_df %>% 
    dplyr::mutate(
      prob_buy_sim = NA,
      demand = NA,
      revenue = NA
    )
  
  # t期における消費者の分布、購買割合、需要、収入を計算
  for (t in 1:num_period_sim) {
    # 初期、かつ、一時的な値下げの場合にのみ、一時的な値下げのCCPを参照する
    if (t == 1 & scenario == 'Temporary'){
      # 分布の推移
      # 注意: 最終期は計算する必要がない
      if (t != num_period_sim){
        consumer_dist_sim[t+1,] <- 
          consumer_dist_sim[t,] %*% G_CCP_list[['Temporary']]
      }
      # 購買割合の計算
      discount_sim_df$prob_buy_sim[t] <- 
        consumer_dist_sim[t,] %*% CCP_list[['Temporary']][,2]
      # 需要量の計算
      discount_sim_df$demand[t] <-
        discount_sim_df$prob_buy_sim[t] * num_consumer_sim
      # 収入の計算 
      discount_sim_df$revenue[t] <- 
        state_df %>% 
        dplyr::mutate(consumer_dist = consumer_dist_sim[t,],
                      CCP = CCP_list[['Temporary']][,2],
                      revenue = (price - 100) * consumer_dist * CCP * num_consumer_sim) %>% 
        dplyr::summarise(revenue = sum(revenue)) %>% 
        as.matrix()
    } else{
      # 実際に消費者が従うCCPのシナリオ名を得る
      # 注意: 一時的な値下げの場合、t >= 2においてはベースラインのCCPを参照する
      scenario_current <- 
        case_when(scenario %in% c('Baseline', 'Permanent') ~ scenario,
                  scenario == 'Temporary' ~ 'Baseline')
      # 分布の推移
      # 注意: 最終期は計算する必要がない
      if (t != num_period_sim){
        consumer_dist_sim[t+1,] <- 
          consumer_dist_sim[t,] %*% G_CCP_list[[scenario_current]]
      }
      # 購買割合の計算
      discount_sim_df$prob_buy_sim[t] <- 
        consumer_dist_sim[t,] %*% CCP_list[[scenario_current]][,2]
      # 需要量の計算
      discount_sim_df$demand[t] <-
        discount_sim_df$prob_buy_sim[t] * num_consumer_sim
      # 収入の計算
      # 値下げをしている場合のみ、値下げ額を取得
      discount <- 
        case_when(scenario %in% c('Baseline', 'Temporary') ~ 0,
                  scenario == 'Permanent' ~ 100)
      discount_sim_df$revenue[t] <- 
        state_df %>% 
        dplyr::mutate(consumer_dist = consumer_dist_sim[t,],
                      CCP = CCP_list[[scenario_current]][,2],
                      revenue = (price - discount) * consumer_dist * CCP * num_consumer_sim) %>% 
        dplyr::summarise(revenue = sum(revenue)) %>% 
        as.matrix()
    }
  }
  # 購買割合、需要、収入の変数名をシナリオに応じて変更
  discount_sim_df <- 
    discount_sim_df %>% 
    dplyr::rename(!!str_c('prob_buy_sim_', scenario) := 'prob_buy_sim',
                  !!str_c('demand_', scenario) := 'demand',
                  !!str_c('revenue_', scenario) := 'revenue')
}

7.3 値下げの効果

# ベースラインを基準とした、購買割合の変化を計算
discount_sim_df %>%
  dplyr::mutate(PermanentChange = prob_buy_sim_Permanent - prob_buy_sim_Baseline,
                TemporaryChange = prob_buy_sim_Temporary - prob_buy_sim_Baseline) %>%
  dplyr::select(period, PermanentChange, TemporaryChange) %>%
  tidyr::pivot_longer(cols = -period) %>%
  dplyr::rename(scenario = name, prob_buy = value) %>%
  ggplot2::ggplot(aes(x = period, y = prob_buy, color = scenario)) +
  geom_line() +
  geom_point() +
  ggplot2::theme_minimal() +
  # ggplot2::scale_x_continuous(breaks = seq(0, num_period_sim, by = 1)) +
  ggplot2::geom_hline(yintercept = 0,
                      linetype = 'dashed',
                      color = "blue", size = 0.3)

8 需要と収入の計算

今回の反実仮想で扱った全てのシナリオ（EDLPと値下げ）の需要と収入を計算する。そのために、まずEDLPの場合の消費者の分布をシミュレートする。

8.1 消費者の分布のシミュレーション（EDLP）

# データから状態変数毎の消費者の分布を計算。
# mileage_idの順番に並んだ、(1 * 21) の行ベクトル。
consumer_dist_obs_edlp <- 
  data_gen %>% 
  dplyr::group_by(mileage_id) %>% 
  dplyr::summarise(num_obs = n(),
                   .groups = 'drop') %>% 
  dplyr::mutate(consumer_dist_obs_edlp = num_obs/sum(num_obs)) %>% 
  dplyr::arrange(mileage_id) %>% 
  dplyr::select(consumer_dist_obs_edlp) %>% 
  as.matrix() %>% 
  t()

# EDLPのシミュレーション結果を表すデータフレームを前もって作成
edlp_sim_df <-
  dplyr::tibble(period = 1:num_period_sim)

# すべての価格についてシミュレーションを行う
for (fixed_price in price_states) {
  scenario_edlp <- stringr::str_c('edlp', fixed_price)
  
  # 消費者の分布の推移を表す行列を前もって作成
  consumer_dist_sim <- matrix(0, num_period_sim, length(consumer_dist_obs_edlp))
  # 初期の分布をデータ上の分布とする
  consumer_dist_sim[1,] <- consumer_dist_obs_edlp
  
  # 購買割合、需要、収入の空の列を作成
  edlp_sim_df <-
    edlp_sim_df %>% 
    dplyr::mutate(
      prob_buy_sim = NA,
      demand = NA,
      revenue = NA
    )
  
  # t期における消費者の分布、購買割合、需要、収入を計算
  for (t in 1:num_period_sim) {
    # 分布の推移
      # 注意: 最終期は計算する必要がない
    if (t != num_period_sim){
      consumer_dist_sim[t+1,] <- 
        consumer_dist_sim[t,] %*% G_CCP_list[[scenario_edlp]]
    }
    # 購買割合の計算
    edlp_sim_df$prob_buy_sim[t] <- 
      consumer_dist_sim[t,] %*% CCP_list[[scenario_edlp]][,2]
    # 需要量の計算
    edlp_sim_df$demand[t] <-
      edlp_sim_df$prob_buy_sim[t] * num_consumer_sim
    # 収入の計算
    edlp_sim_df$revenue[t] <- 
      state_df %>% 
      dplyr::filter(price == fixed_price) %>% 
      dplyr::arrange(mileage_id) %>% 
      dplyr::mutate(consumer_dist = consumer_dist_sim[t,],
                    CCP = CCP_list[[scenario_edlp]][,2],
                    revenue = fixed_price * consumer_dist * CCP * num_consumer_sim) %>% 
      dplyr::summarise(revenue = sum(revenue)) %>% 
      as.matrix()
  }
  # 購買割合、需要、収入の変数名をシナリオに応じて変更
  edlp_sim_df <- 
    edlp_sim_df %>% 
    dplyr::rename(!!str_c('prob_buy_sim_', scenario_edlp) := 'prob_buy_sim',
                  !!str_c('demand_', scenario_edlp) := 'demand',
                  !!str_c('revenue_', scenario_edlp) := 'revenue')
}

# 結果を図示するためにすべてのシミュレーション結果を統合
sim_df <- 
  discount_sim_df %>% 
  dplyr::left_join(edlp_sim_df, by = 'period')

8.2 シナリオ毎の需要と収入の比較

8.2.1 EDLPの需要と収入の推移

sim_df %>% 
  dplyr::select(period, stringr::str_c('demand_edlp', price_states)) %>% 
  tidyr::pivot_longer(cols = - period) %>% 
  dplyr::rename(scenario = name, demand = value) %>% 
  ggplot2::ggplot(aes(x = period, y = demand, color = scenario)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() +
  ggplot2::theme_minimal()

sim_df %>% 
  dplyr::select(period, stringr::str_c('revenue_edlp', price_states)) %>% 
  tidyr::pivot_longer(cols = - period) %>% 
  dplyr::rename(scenario = name, revenue = value) %>% 
  ggplot2::ggplot(aes(x = period, y = revenue, color = scenario)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() +
  ggplot2::theme_minimal()

8.2.2 値下げをした場合の需要と収入の推移

# 需要の推移を図示
sim_df %>% 
  dplyr::select(period, stringr::str_c('demand_', discount_scenerio_vec)) %>% 
  tidyr::pivot_longer(cols = - period) %>% 
  dplyr::rename(scenario = name, demand = value) %>% 
  ggplot2::ggplot(aes(x = period, y = demand, color = scenario)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() + 
  ggplot2::theme_minimal()

# 需要の累積和を図示
sim_df %>% 
  dplyr::filter(period <= 20) %>% 
  dplyr::select(period, stringr::str_c('demand_', discount_scenerio_vec)) %>% 
  dplyr::mutate(across(starts_with('demand'), cumsum)) %>% 
  tidyr::pivot_longer(cols = - period) %>% 
  dplyr::rename(scenario = name, demand = value) %>% 
  ggplot2::ggplot(aes(x = period, y = demand, color = scenario)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() + 
  ggplot2::theme_minimal()

# ベースラインを100とした時の、需要の累積和を図示
sim_df %>% 
  dplyr::filter(period <= 20) %>% 
  dplyr::select(period, stringr::str_c('demand_', discount_scenerio_vec)) %>% 
  dplyr::mutate(across(starts_with('demand'), cumsum)) %>% 
  dplyr::mutate(Permanent = demand_Permanent/demand_Baseline * 100,
                Temporary = demand_Temporary/demand_Baseline * 100) %>%
  dplyr::select(period, Permanent, Temporary) %>%
  tidyr::pivot_longer(cols = - period) %>% 
  dplyr::rename(scenario = name, demand = value) %>% 
  ggplot2::ggplot(aes(x = period, y = demand, color = scenario)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() + 
  ggplot2::theme_minimal()

sim_df %>% 
  dplyr::select(period, stringr::str_c('revenue_', discount_scenerio_vec)) %>% 
  tidyr::pivot_longer(cols = - period) %>% 
  dplyr::rename(scenario = name, revenue = value) %>% 
  ggplot2::ggplot(aes(x = period, y = revenue, color = scenario)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() +
  ggplot2::theme_minimal()

8.2.3 合計の需要と収入

demand_df <-
  sim_df %>% 
  dplyr::summarise(across(starts_with('demand'), sum)) %>% 
  dplyr::rename_with(~str_replace(., "demand_", "")) %>% 
  tidyr::pivot_longer(cols = everything()) %>% 
  dplyr::rename(scenario = name, demand = value)

revenue_df <- 
  sim_df %>% 
  # 割引現在価値を計算するため、割引因子の列を追加
  dplyr::mutate(beta = beta ^ (period-1)) %>% 
  dplyr::summarise(across(starts_with('revenue'), ~sum(.x * beta))) %>% 
  dplyr::rename_with(~str_replace(., "revenue_", "")) %>% 
  tidyr::pivot_longer(cols = everything()) %>% 
  dplyr::rename(scenario = name, revenue = value)

demand_df %>% 
  dplyr::left_join(revenue_df, by = 'scenario')%>% 
  dplyr::mutate( per_rev = revenue/demand/10) %>% 
  kableExtra::kable(format = "html", 
                    booktabs = TRUE) %>% 
  kableExtra::kable_styling()

scenario	demand	revenue	per_rev
Baseline	579.1927	1190089	205.4738
Permanent	620.0550	1219207	196.6288
Temporary	582.0691	1193790	205.0941
edlp2000	709.7349	1296483	182.6714
edlp2100	661.4560	1267340	191.5985
edlp2200	617.3481	1237757	200.4958
edlp2300	576.6884	1207384	209.3650
edlp2400	538.9096	1175943	218.2079
edlp2500	503.5640	1143223	227.0263

経済セミナー連載「実証ビジネス・エコノミクス」
第８回「価格戦略をダイナミックに考える：シングルエージェント動学モデルの推定［応用編］」

上武康亮・遠山祐太・若森直樹・渡辺安虎

最終更新: 2022-05-24

1 はじめに

1.1 留意事項

1.2 謝辞

1.3 全体の流れ

2 下準備

2.1 Rに関する下準備

2.2 データの読み込み

3 記述統計の確認

4 二段階推定方法の実践

4.1 Step 1: CCP の推定、及び State transition の推定

4.1.1 State transition の推定

4.1.2 CCPの推定

4.2 Step 2: パラメタの推定

4.2.1 行列形式によるインバージョン

4.2.2 有限依存性 (finite dependence) アプローチ

5 推定方法の比較

5.1 第7回の復習：入れ子不動点アルゴリズム

5.2 推定方法の比較表

6 反実仮想分析1: Every Day Low Price

6.1 EDLPの効果

7 反実仮想分析2: 永続的・一時的な値下げ

7.1 値下げが起きた時のCCPの計算

7.2 消費者の分布のシミュレーション（値下げ）

7.3 値下げの効果

8 需要と収入の計算

8.1 消費者の分布のシミュレーション（EDLP）

8.2 シナリオ毎の需要と収入の比較

8.2.1 EDLPの需要と収入の推移

8.2.2 値下げをした場合の需要と収入の推移

8.2.3 合計の需要と収入

経済セミナー 連載「実証ビジネス・エコノミクス」 第８回 「価格戦略をダイナミックに考える：シングルエージェント動学モデルの推定［応用編］」

上武康亮・遠山祐太・若森直樹・渡辺安虎

最終更新: 2022-05-24

1 はじめに

1.1 留意事項

1.2 謝辞

1.3 全体の流れ

2 下準備

2.1 Rに関する下準備

2.2 データの読み込み

3 記述統計の確認

4 二段階推定方法の実践

4.1 Step 1: CCP の推定、及び State transition の推定

4.1.1 State transition の推定

4.1.2 CCPの推定

4.2 Step 2: パラメタの推定

4.2.1 行列形式によるインバージョン

4.2.2 有限依存性 (finite dependence) アプローチ

5 推定方法の比較

5.1 第7回の復習：入れ子不動点アルゴリズム

5.2 推定方法の比較表

6 反実仮想分析1: Every Day Low Price

6.1 EDLPの効果

7 反実仮想分析2: 永続的・一時的な値下げ

7.1 値下げが起きた時のCCPの計算

7.2 消費者の分布のシミュレーション（値下げ）

7.3 値下げの効果

8 需要と収入の計算

8.1 消費者の分布のシミュレーション（EDLP）

8.2 シナリオ毎の需要と収入の比較

8.2.1 EDLPの需要と収入の推移

8.2.2 値下げをした場合の需要と収入の推移

8.2.3 合計の需要と収入

経済セミナー連載「実証ビジネス・エコノミクス」
第８回「価格戦略をダイナミックに考える：シングルエージェント動学モデルの推定［応用編］」