経済セミナー 連載「実証ビジネス・エコノミクス」
第７回 「将来予想のインパクトを測る：シングルエージェント動学モデルの推定［基礎編］」

1 はじめに

2022年4月・5月号の「実証ビジネス・エコノミクス」の第7回「将来予想のインパクトを測る：シングルエージェント動学モデルの推定［基礎編］」に付随するRコードの紹介になります。

1.1 留意事項

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

1.2 謝辞

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

1.3 全体の流れ

1. 下準備：パッケージの導入、パラメータの設定、データの生成
2. 記述統計の作成
3. 遷移行列の推定
4. Nested Fixed Point Algorithmによるパラメータの推定

2 下準備：パッケージの導入、パラメータの設定、データの生成

2.1 Rに関する下準備

今回はevd,numDeriv,plot3Dを新たに利用する。 なお、以下のrequireで読み込んでいるパッケージについてまだインストールしていない場合には、最初にインストールすること。

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

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

## Warning: パッケージ 'skimr' はバージョン 4.1.2 の R の下で造られました

# install.packages("numDeriv")
require(evd)

## Warning: パッケージ 'evd' はバージョン 4.1.3 の R の下で造られました

# install.packages("evd")
require(numDeriv)
# install.packages("plot3D")
library("plot3D")

## Warning: パッケージ 'plot3D' はバージョン 4.1.3 の R の下で造られました

2.2 データの生成

今回は現実のデータではなく、仮想データを使用する。そこで真のパラメータや遷移行列を設定し、データを生成していく。

2.2.1 パラメータの設定

## パラメータの設定

# 走行距離にかかるtheta_cと車の価格にかかるtheta_p
theta_true <- c(theta_c = 0.004, theta_p = 0.003)

# 時間割引率
beta <- 0.99

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

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

2.2.1.1 状態変数の作成

## 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 x 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

2.2.1.2 遷移行列の作成

今回、走行距離と価格が遷移確率によって変化していく。そこでそれぞれの遷移確率から遷移行列を作成する。

走行距離の遷移行列を作成する関数

# パラメタを所与として走行距離の遷移行列を出力する関数を作成
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)
}

走行距離の遷移行列を作成

# 走行距離の遷移行列のパラメタを設定し、遷移行列を作成する
kappa_true <- c(0.25, 0.05)

mileage_trans_mat_true <- gen_mileage_trans(kappa_true)

# 走行距離の遷移行列の4行4列までを表示
mileage_trans_mat_true[1:4,1:4,1]

##      [,1] [,2] [,3] [,4]
## [1,]  0.7 0.25 0.05 0.00
## [2,]  0.0 0.70 0.25 0.05
## [3,]  0.0 0.00 0.70 0.25
## [4,]  0.0 0.00 0.00 0.70

価格の遷移行列を作成

# 価格の遷移行列のパラメタを設定し、遷移行列を作成する
lambda_true <- c(0.1, 0.2, 0.2, 0.2, 0.2,
                 0.1, 0.2, 0.2, 0.2, 0.2,
                 0.1, 0.1, 0.2, 0.2, 0.1,
                 0.1, 0.1, 0.2, 0.2, 0.1,
                 0.05, 0.05, 0.1, 0.1, 0.2,
                 0.05, 0.05, 0.1, 0.1, 0.2)

price_trans_mat_true <- gen_price_trans(lambda_true)

# 価格の遷移行列を表示
price_trans_mat_true

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.10 0.10  0.2  0.2  0.2  0.2
## [2,] 0.10 0.10  0.2  0.2  0.2  0.2
## [3,] 0.10 0.10  0.3  0.2  0.2  0.1
## [4,] 0.10 0.10  0.2  0.3  0.2  0.1
## [5,] 0.05 0.05  0.1  0.1  0.5  0.2
## [6,] 0.05 0.05  0.1  0.1  0.2  0.5

# コントロール変数毎の遷移行列を作成
trans_mat_true <- list()

# 車を購入しない場合の遷移行列
trans_mat_true$not_buy <- 
  mileage_trans_mat_true[,,1] %x% price_trans_mat_true

# 車を購入する場合の遷移行列
trans_mat_true$buy <- 
  mileage_trans_mat_true[,,2] %x% price_trans_mat_true

# 定常状態での価格の分布を計算
# 以下を満たすような price_dist_steady を求める
# price_dist_steady %*% price_trans_mat == price_dist_steady

# 固有値/固有ベクトルを求める
# 固有値が1となる固有ベクトルは1つだけ（1つめ）
price_trans_eigen <- eigen(t(price_trans_mat_true))

# 価格の定常分布を求める
price_dist_steady <-
  price_trans_eigen$vectors[,1]/sum(price_trans_eigen$vectors[,1])

price_dist_steady

## [1] 0.07377049 0.07377049 0.16393443 0.16393443 0.28571429 0.23887588

2.2.1.3 効用関数の定義

効用関数は以下の式で定義される。 \[ u(x_t, i_t; \theta) = \begin{cases} - \theta_c m_t & \text{if} \quad i_t = 0 \\ - \theta_p p_t & \text{if} \quad i_t = 1 \end{cases} \]

# 状態変数、コントロール変数毎の今期の効用を返す関数
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)
}

2.2.1.4 価値関数反復法

contraction <- 
  function(theta, beta, trans_mat, state_df) {
    # パラメタより今期の効用を計算
    U <- flow_utility(theta, state_df)
    
    # 価値関数の初期値
    EV_old <- matrix(0, nrow = num_states, ncol = num_choice)
    
    # 価値関数の差の初期値
    diff <- 1000
    
    # 縮小写像の誤差範囲
    tol_level <- 1.0e-10
    
    while (diff > tol_level) {
      # 選択ごとの価値関数を計算
      EV_new <- cbind(
        EV_not_buy <- 
          Euler_const + trans_mat$not_buy %*% log(rowSums(exp(U + beta*EV_old))),
        EV_buy <-
          Euler_const + trans_mat$buy %*% log(rowSums(exp(U + beta*EV_old)))
      )
      # 価値関数の更新による差を評価
      diff <- sum(abs(EV_new-EV_old))
      
      # 価値関数を次のループに渡す(価値関数の更新)
      EV_old <- EV_new
    }
    EV <- EV_old
    colnames(EV) <- c("EV_not_buy", "EV_buy")
    return(EV)
  }

# EVを求める
start_time <- proc.time()

EV_true <- contraction(theta_true, beta, trans_mat_true, state_df)

end_time <- proc.time()
cat("Runtime:\n")

## Runtime:

print((end_time - start_time)[[3]])

## [1] 0.15

# 選択毎の価値関数を定義する
U_true <- flow_utility(theta_true, state_df)
V_CS_true <- U_true + beta*EV_true
colnames(V_CS_true) <- c("V_not_buy", "V_buy")

# state(p,m)ごとに、logitで計算される理論上の条件付き購入確率を計算
prob_buy_true_mat <- matrix(exp(V_CS_true[,"V_buy"])/rowSums(exp(V_CS_true)), 
                            nrow = num_price_states, ncol = num_mileage_states)
prob_buy_true_mat

##              [,1]         [,2]        [,3]        [,4]        [,5]        [,6]
## [1,] 0.0024726232 0.0042603618 0.007049664 0.011200433 0.017090266 0.025063775
## [2,] 0.0018329389 0.0031596426 0.005232079 0.008321642 0.012717111 0.018689107
## [3,] 0.0013585200 0.0023423228 0.003879975 0.006174285 0.009442589 0.013891293
## [4,] 0.0010067708 0.0017363494 0.002877500 0.004582101 0.007014286 0.010332322
## [5,] 0.0007460288 0.0012874740 0.002135678 0.003405593 0.005223411 0.007714133
## [6,] 0.0005527786 0.0009541767 0.001583334 0.002526076 0.003877171 0.005731449
##             [,7]       [,8]       [,9]      [,10]      [,11]      [,12]
## [1,] 0.035376761 0.04815055 0.06334970 0.08078789 0.10015701 0.12106774
## [2,] 0.026450272 0.03612160 0.04771404 0.06112911 0.07617557 0.09259467
## [3,] 0.019687178 0.02693284 0.03565272 0.04579239 0.05722911 0.06978856
## [4,] 0.014667999 0.02010868 0.02668660 0.03437726 0.04310577 0.05275758
## [5,] 0.010987109 0.01512297 0.02016506 0.02611673 0.03294395 0.04058092
## [6,] 0.008173405 0.01126772 0.01505291 0.01953904 0.02470914 0.03052275
##           [,13]      [,14]      [,15]      [,16]      [,17]      [,18]
## [1,] 0.14308885 0.16577560 0.18868129 0.21134866 0.23327754 0.25385965
## [2,] 0.11008546 0.12832312 0.14696556 0.16564457 0.18393734 0.20130915
## [3,] 0.08326164 0.09741644 0.11200146 0.12673711 0.14129141 0.15523193
## [4,] 0.06318930 0.07423618 0.08571309 0.09740637 0.10905285 0.12029899
## [5,] 0.04893600 0.05789523 0.06732107 0.07704365 0.08684085 0.09640043
## [6,] 0.03691954 0.04382114 0.05112915 0.05871736 0.06641473 0.07397364
##           [,19]     [,20]      [,21]
## [1,] 0.27227971 0.2872751 0.29731191
## [2,] 0.21702523 0.2299391 0.23864342
## [3,] 0.16795126 0.1784896 0.18565187
## [4,] 0.13063724 0.1392574 0.14514098
## [5,] 0.10526151 0.1126835 0.11773586
## [6,] 0.08102088 0.0869503 0.09099472

2.2.2 シミュレーション

今回は1000人の消費者が50年にわたって車を購入するか、しないかの意思決定をしてきたことを想定したデータを作成する。ただし、私たちが観測できるのは最後の10年とする。

## サンプルサイズを決める

# 1000人の消費者が存在
num_consumer <- 1000

# 50年分の月次データを生成した後、最後の10年のみが観察できるとする
num_period <- 12 * 50
num_period_obs <- 12 * 10

# 総観察数
num_obs <- num_consumer * num_period

# 累積分布確率を持つように遷移行列を変換（行方向に足し上げる）
trans_mat_cum <- list()
trans_mat_cum$not_buy <- t(apply(trans_mat_true$not_buy, 1, cumsum))
trans_mat_cum$buy <- t(apply(trans_mat_true$buy, 1, cumsum))

# 乱数を固定
set.seed(1)

# 生成するデータの元となるdata.frameを作成
data_gen <- 
  dplyr::tibble(
    consumer = rep(1:num_consumer, each = num_period),
    period = rep(1:num_period, times = num_consumer),
    eps_type1_not_buy = evd::rgev(num_obs),
    eps_type1_buy = evd::rgev(num_obs),
    eps_unif = runif(num_obs),
    eps_price_state_unif = runif(num_obs),
    state_id = 0,
    action = 0
  )

# 各消費者についてデータを生成する関数を作成
generate_data <- function(df, V_CS, state_df, price_dist_steady) {
  
  # Step 1: 各消費者について、初期のstate_idを決める
  # 価格は定常分布に従うとし、走行距離は0とする
  
  # 価格の定常分布の累積値を計算
  price_dist_steady_cumsum <- cumsum(price_dist_steady)
  
  # 一様分布から生成した乱数が、定常分布の累積値を
  # 初めて下回ったところを1期の状態変数とする
  price_id_consumer <- 0
  exceed_trans_prob_price <- TRUE
  while(exceed_trans_prob_price) {
      price_id_consumer <- price_id_consumer + 1
      exceed_trans_prob_price <- 
        (df$eps_price_state_unif[1] >
           price_dist_steady_cumsum[price_id_consumer])
  }
  
  # state_idに変換し、各消費者の初期のstate_idを決める
  df$state_id[1] <-  state_df %>% 
    # mileageは0とする
    dplyr::filter(mileage_id == 1) %>% 
    dplyr::filter(price_id == price_id_consumer) %>% 
    dplyr::select(state_id) %>% 
    as.numeric()
  
  # Step 2: 各消費者について、状態変数、コントロール変数を逐次的に生成
  for (t in 1:(num_period-1)) {
    # t期のstateを取得
    state_id_today <- df$state_id[t]
    
    # 価値関数に基づいて、購入するかどうかを決める
    if (V_CS[,'V_not_buy'][state_id_today] + df$eps_type1_not_buy[t] > 
        V_CS[,'V_buy'][state_id_today] + df$eps_type1_buy[t]){
      
      # 購入しない
      df$action[t] <- 0
      
      # 直面する遷移行列を定義
      trans_mat_cum_today <- trans_mat_cum$not_buy
      
    }else{
      # 購入する
      df$action[t] <- 1
      
      # 直面する遷移行列を定義
      trans_mat_cum_today <- trans_mat_cum$buy
      
    }
    
    # t+1期のstateを決める
    state_id_tomorrow <- 0
    exceed_trans_prob <- TRUE
    # 一様分布から生成した乱数が、遷移確率の累積分布の値を
    # 初めて下回ったところをt+1期の状態変数とする
    while (exceed_trans_prob) {
      state_id_tomorrow <- state_id_tomorrow + 1
      trans_prob <- trans_mat_cum_today[state_id_today, state_id_tomorrow]
      exceed_trans_prob <- (df$eps_unif[t] > trans_prob)
    }
    df$state_id[t+1]<- state_id_tomorrow
  }
  return(df)
}

data_gen <- 
  data_gen %>%
  # 消費者ごとにデータを分割
  dplyr::group_split(consumer) %>%
  # 上記の関数で定義したデータの生成過程をすべての消費者に対して行う
  purrr::map_dfr(generate_data,
                 V_CS = V_CS_true,
                 state_df = state_df, 
                 price_dist_steady = price_dist_steady) %>% 
  # 最後の10年のみが観察できるとする
  dplyr::filter(period > (num_period - num_period_obs)) %>% 
  # 状態変数を表した列を追加
  dplyr::left_join(state_df, by = 'state_id')

data_gen %>% tail(3)

## # A tibble: 3 x 12
##   consumer period eps_type1_not_buy eps_type1_buy eps_unif eps_price_state_unif
##      <int>  <int>             <dbl>         <dbl>    <dbl>                <dbl>
## 1     1000    598            -1.08          0.264   0.212                 0.239
## 2     1000    599             0.639         0.501   0.0818                0.348
## 3     1000    600             0.776         4.62    0.235                 0.969
## # ... with 6 more variables: state_id <dbl>, action <dbl>, price_id <int>,
## #   mileage_id <int>, price <dbl>, mileage <dbl>

# 不要なオブジェクトの削除
rm(V_CS_true, trans_mat_cum)

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

走行距離と価格の分布を確認する

# 走行距離の分布
data_gen %>%
  ggplot(aes(x = price)) + geom_histogram(binwidth = 100)

# 価格の分布
data_gen %>%
  ggplot(aes(x = mileage)) + geom_histogram(binwidth = 5)

走行距離と購入パターンの関連

# それぞれの走行距離において購入した割合を観察
data_gen %>% 
  dplyr::group_by(mileage) %>% 
  dplyr::summarize(num_state = n(),
                   sum_action = sum(action)) %>% 
  dplyr::mutate(prob_buy = sum_action / num_state) %>% 
  ggplot(aes(x = mileage, y = prob_buy)) + 
  geom_bar(stat = "identity")

価格と購入パターンの関連

# それぞれの価格において購入した割合を観察
data_gen %>% 
  dplyr::group_by(price) %>% 
  dplyr::summarize(num_state = n(),
                   sum_action = sum(action),
                   .groups = 'drop') %>% 
  dplyr::mutate(prob_buy = sum_action / num_state) %>% 
  ggplot(aes(x = price, y = prob_buy)) + 
  geom_bar(stat = "identity")

観測された条件付き購入確率

# state(p,m)ごとに、観測された条件付き購入確率を計算
prob_buy_obs_mat <- 
  data_gen %>%
  dplyr::group_by(mileage,price) %>%
  dplyr::summarize(num_state = n(),
                   sum_action = sum(action),
                   .groups = 'drop') %>%
  dplyr::mutate(prob_buy = sum_action / num_state) %>% 
  dplyr::select(prob_buy) %>% 
  as.matrix() %>% 
  matrix(nrow = num_price_states, ncol = num_mileage_states)
prob_buy_obs_mat

##              [,1]         [,2]        [,3]        [,4]        [,5]        [,6]
## [1,] 0.0051107325 0.0066755674 0.002724796 0.005952381 0.014814815 0.020527859
## [2,] 0.0000000000 0.0013089005 0.003680982 0.014965986 0.008902077 0.021613833
## [3,] 0.0007390983 0.0019379845 0.005660377 0.005901639 0.009870450 0.011118833
## [4,] 0.0007423905 0.0006309148 0.001244555 0.003076923 0.005685407 0.014512785
## [5,] 0.0004621072 0.0007256894 0.002724796 0.003183587 0.004368402 0.006627393
## [6,] 0.0000000000 0.0017969452 0.001351960 0.002937474 0.004672897 0.003872633
##            [,7]       [,8]       [,9]      [,10]      [,11]      [,12]
## [1,] 0.02579666 0.03826955 0.07155635 0.07441016 0.10675381 0.13065327
## [2,] 0.02567976 0.03271028 0.04887218 0.05313093 0.08597285 0.08418367
## [3,] 0.01908657 0.02688953 0.03476821 0.04000000 0.06701571 0.07069555
## [4,] 0.01562500 0.02339181 0.02795031 0.03206751 0.04459203 0.06221719
## [5,] 0.01140251 0.01518260 0.02557201 0.02721425 0.03360000 0.03564223
## [6,] 0.01204819 0.01099450 0.02508179 0.02009185 0.02197071 0.02961808
##           [,13]      [,14]      [,15]      [,16]      [,17]      [,18]
## [1,] 0.11235955 0.14482759 0.16017316 0.20105820 0.26271186 0.24242424
## [2,] 0.09969789 0.15267176 0.13425926 0.13609467 0.17037037 0.17391304
## [3,] 0.06250000 0.08239095 0.12173913 0.11750600 0.12786885 0.21212121
## [4,] 0.07061791 0.06060606 0.08823529 0.10169492 0.10144928 0.11518325
## [5,] 0.03657588 0.06717124 0.06616052 0.05898876 0.08266129 0.10344828
## [6,] 0.03083333 0.04535147 0.06342780 0.05972696 0.05239180 0.06188925
##           [,19]      [,20]      [,21]
## [1,] 0.25396825 0.24528302 0.22222222
## [2,] 0.19230769 0.31666667 0.17355372
## [3,] 0.17361111 0.15596330 0.17727273
## [4,] 0.13669065 0.08695652 0.17427386
## [5,] 0.09722222 0.11797753 0.12831858
## [6,] 0.08298755 0.03750000 0.09917355

hist3D(x = mileage_states, y = price_states, z = t(prob_buy_obs_mat), zlim=c(0,0.4),
       bty = "g", phi = 10,  theta = -60, axes=TRUE,label=TRUE,
        xlab = "Mileage", ylab = "Price", zlab = "Probability", main = "Conditional probability of buying",
        col = "#0080ff", border = "blue", shade = 0.4,
        ticktype = "detailed", space = 0.05, d = 2, cex.axis = 0.8)

4 遷移行列の推定

ここでは遷移行列をデータから推定する

4.1 走行距離の遷移行列の推定

# 遷移行列の推定で使うため、ラグ変数を追加
data_gen <- 
  data_gen %>% 
  dplyr::group_by(consumer) %>% 
  dplyr::mutate(lag_price_id = lag(price_id),
                lag_mileage_id = lag(mileage_id),
                lag_action = lag(action)) %>% 
  dplyr::ungroup() 

# それぞれの確率が実現した観察の数を数える
num_cond_obs_mileage <- 
  data_gen %>% 
  # 1期目は推定に使えないため落とす
  dplyr::filter(period != (num_period - num_period_obs + 1)) %>% 
  # 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]))

# 最尤法の解析解により標準誤差を求める
Infomat_mileage_est <- matrix(0, nrow = 2, ncol = 2)

# 最尤法のフィッシャー情報量を求める
Infomat_mileage_est[1,1] <- 
  (num_cond_obs_mileage[1] / (1 - kappa_est[1] - kappa_est[2])^2) +
    (num_cond_obs_mileage[2] / kappa_est[1]^2) +
    (num_cond_obs_mileage[4] / (kappa_est[1]+kappa_est[2])^2)
Infomat_mileage_est[1,2] <- 
  (num_cond_obs_mileage[1] / (1 - kappa_est[1] - kappa_est[2])^2) +
    (num_cond_obs_mileage[4] / (kappa_est[1]+kappa_est[2])^2)
Infomat_mileage_est[2,1] <- Infomat_mileage_est[1,2]
Infomat_mileage_est[2,2] <- 
  (num_cond_obs_mileage[1] / (1 - kappa_est[1] - kappa_est[2])^2) +
    (num_cond_obs_mileage[3] / kappa_est[2]^2) +
    (num_cond_obs_mileage[4] / (kappa_est[1]+kappa_est[2])^2)

# 逆行列の対角要素の平方根が標準誤差になる
kappa_se <- sqrt(diag(solve(Infomat_mileage_est)))

dplyr::tibble(kappa_est, kappa_se)

## # A tibble: 2 x 2
##   kappa_est kappa_se
##       <dbl>    <dbl>
## 1    0.251  0.00126 
## 2    0.0500 0.000637

4.2 価格の遷移行列の推定

# それぞれの確率が実現した観察の数を数える
num_cond_obs_price <- 
  data_gen %>% 
  # 1期目は推定に使えないため落とす
  dplyr::filter(period != (num_period - num_period_obs + 1)) %>% 
  # 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_se <- c()
for (i in 1:num_price_states) {
  # 最尤法のフィッシャー情報量を求める
  Infomat_price_est <- 
    diag(num_cond_obs_price[i,],
         num_price_states)[-i,-i] / 
    (lambda_est_mat[-i,-i] ^ 2) + 
    (num_cond_obs_price[i,i] / 
       lambda_est_mat[i,i] ^ 2) *
    matrix(1, num_price_states, num_price_states)[-i,-i]
  lambda_se <- c(
    lambda_se,
    # 逆行列の対角要素の平方根が標準誤差になる
    sqrt(diag(solve(Infomat_price_est)))
  )
}

lambda_se_mat <- 
    c(0, lambda_se[1], lambda_se[2], lambda_se[3], lambda_se[4], lambda_se[5],
      lambda_se[6], 0, lambda_se[7], lambda_se[8], lambda_se[9], lambda_se[10],
      lambda_se[11], lambda_se[12], 0, lambda_se[13], lambda_se[14], lambda_se[15],
      lambda_se[16], lambda_se[17], lambda_se[18], 0, lambda_se[19], lambda_se[20],
      lambda_se[21], lambda_se[22], lambda_se[23], lambda_se[24], 0, lambda_se[25],
      lambda_se[26], lambda_se[27], lambda_se[28], lambda_se[29], lambda_se[30], 0) %>% 
    matrix(ncol = num_price_states, nrow = num_price_states, byrow=T)
lambda_se_mat

##             [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
## [1,] 0.000000000 0.003244047 0.006662820 0.006789639 0.009766072 0.009609858
## [2,] 0.003391513 0.000000000 0.006731948 0.006603122 0.009625120 0.009594420
## [3,] 0.002303224 0.002236035 0.000000000 0.004685278 0.006898800 0.007832604
## [4,] 0.002336093 0.002220141 0.004597027 0.000000000 0.006800609 0.007718928
## [5,] 0.002430599 0.002334927 0.004489448 0.004562071 0.000000000 0.004960462
## [6,] 0.002633243 0.002530651 0.004935344 0.004988288 0.005458357 0.000000000

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

## # A tibble: 30 x 2
##    lambda_est lambda_se
##         <dbl>     <dbl>
##  1      0.104   0.00324
##  2      0.199   0.00666
##  3      0.199   0.00679
##  4      0.192   0.00977
##  5      0.203   0.00961
##  6      0.101   0.00339
##  7      0.192   0.00673
##  8      0.210   0.00660
##  9      0.199   0.00963
## 10      0.199   0.00959
## # ... with 20 more rows

5 パラメータの推定

5.1 静学的なロジットによる推定

ここでは消費者が動学的にではなく静学的に意思決定を行っていると仮定し、単純なロジットモデルを用いて\(\theta_c\)と\(\theta_p\)を推定する。

行列の操作に必要な関数を準備する。

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

推定に用いる目的関数を用意する。

# 対数尤度関数を定義
logLH_stat <- function(theta, state_df, df){
  
  
  # 選択毎の効用関数を求める
  U <- flow_utility(theta, state_df)
  # 選択確率を計算
  prob_C_stat <- exp(U) / rowSums(exp(U))
  # 対数尤度を計算
  sum(log(mat_ij(df$state_id, df$action + 1, prob_C_stat)))
}

パラメータ\(\theta_c\)と\(\theta_p\)を推定する

start_time <- proc.time()

# 最適化
logit_stat_opt <- optim(theta_true, logLH_stat,
                  state_df = state_df, df = data_gen, 
                  control = list(fnscale = -1), 
                  method = "Nelder-Mead")

end_time <- proc.time()
cat("Runtime:\n")

## Runtime:

print((end_time - start_time)[[3]])

## [1] 35.7

theta_est_stat <- logit_stat_opt$par
theta_est_stat

##     theta_c     theta_p 
## 0.042051810 0.002368664

標準誤差を推定する

hessian_stat <- numDeriv::hessian(func = logLH_stat, x = theta_est_stat, 
                             state_df = state_df, df = data_gen)
theta_se_stat <- sqrt(diag(solve(-hessian_stat)))
dplyr::tibble(theta_est_stat, theta_se_stat)

## # A tibble: 2 x 2
##   theta_est_stat theta_se_stat
##            <dbl>         <dbl>
## 1        0.0421      0.000687 
## 2        0.00237     0.0000191

5.2 不動点アルゴリズムによる推定

ここでは不動点アルゴリズムを用いて\(\theta_c\)と\(\theta_p\)を推定する。

# 推定された遷移行列を取得
trans_mat_hat <- list()
trans_mat_hat$not_buy <- 
  gen_mileage_trans(kappa_est)[,,1] %x% gen_price_trans(lambda_est)
trans_mat_hat$buy <- 
  gen_mileage_trans(kappa_est)[,,2] %x% gen_price_trans(lambda_est)

推定に用いる目的関数を定義する

# 対数尤度関数を定義
logLH <- function(theta, beta, trans_mat, state_df, df){
  
  # 選択ごとの期待価値関数を計算
  EV <- contraction(theta, beta, trans_mat, state_df)
 
  # 選択毎の価値関数を定義する
  U <- flow_utility(theta, state_df)
  V_CS <- U + beta*EV
  # 選択確率を計算
  prob_C <- exp(V_CS) / rowSums(exp(V_CS))
  # 対数尤度を計算
  sum(log(mat_ij(df$state_id, df$action + 1, prob_C)))
}

パラメータ\(\theta_c\)と\(\theta_p\)を推定する

start_time <- proc.time()

# 最適化
NFXP_opt <- optim(theta_true, logLH,
                  beta = beta, trans_mat = trans_mat_hat, state_df = state_df, df = data_gen, 
                  control = list(fnscale = -1), 
                  method = "Nelder-Mead")

end_time <- proc.time()
cat("Runtime:\n")

## Runtime:

print((end_time - start_time)[[3]])

## [1] 22.72

theta_est <- NFXP_opt$par
theta_est

##     theta_c     theta_p 
## 0.003887129 0.002971684

標準誤差を推定する

hessian <- numDeriv::hessian(func = logLH, x = theta_est, 
                             beta = beta,　trans_mat = trans_mat_hat, state_df = state_df, df = data_gen)
theta_se <- sqrt(diag(solve(-hessian)))
dplyr::tibble(theta_est, theta_se)

## # A tibble: 2 x 2
##   theta_est  theta_se
##       <dbl>     <dbl>
## 1   0.00389 0.0000865
## 2   0.00297 0.0000362