Experimetrics: 17.2 A Level-k Model for the Beauty Contest Game

今回は、実験経済学ではよく知られた美人投票ゲーム(Beauty Contest Game)の架空の被験者データにLevel-kモデルを当てはめ、推論の深さの程度ごとに被験者の割合を推定します。


美人投票ゲームの概要はこちら。

  • 各被験者は0から100までの整数を1つ選ぶ
  • 全被験者が選んだ数字の平均値×(2/3)に最も近い数字を選んだ被験者が勝者

繰り返し支配戦略を消していくことで、最終的には「各被験者が0を選ぶ」というこのゲームのナッシュ均衡にたどり着きます。最初の実験であるNagel (1995)では平均値が35、また多くの被験者が33や22といった数字を選びました。


次にLevel-kモデルの概要です。

  • 人々は様々なレベルの推論を行うと仮定する
  • Level-0(0次)の推論を行うプレイヤー:0から100までの数字を一様にランダムに選択するプレイヤー
  • Level-1(1次)の推論を行うプレイヤー:他のプレイヤーは全てLevel-0の推論を行うとの信念を持ち(つまり、自分以外のプレイヤーが選ぶ数字の平均値は50付近であるとの信念を持つ)、その結果33を最適な推論として選択するプレイヤー
  • Level-2(2次)の推論を行うプレイヤー:他のプレイヤーは全てLevel-1の推論を行うとの信念を持ち(つまり、自分以外のプレイヤーが選ぶ数字の平均値は33付近であるとの信念を持つ)、その結果22を最適な推論として選択するプレイヤー
  • Level-3(3次)の推論を行うプレイヤー:他のプレイヤーは全てLevel-2の推論を行うとの信念を持ち(つまり、自分以外のプレイヤーが選ぶ数字の平均値は22付近であるとの信念を持つ)、その結果15を最適な推論として選択するプレイヤー
  • Level-4(4次)の推論を行うプレイヤー:他のプレイヤーは全てLevel-3の推論を行うとの信念を持ち(つまり、自分以外のプレイヤーが選ぶ数字の平均値は15付近であるとの信念を持つ)、その結果10を最適な推論として選択するプレイヤー
  • 同じ要領でより高次の推論を行うプレイヤーを想定することができるが、今回はLevel-5を最高次とし、他のプレイヤーは全て数字0を選ぶとするナイーブなNashタイプのプレイヤーとする

このように、Level-kモデルによれば、推論がより高次のプレイヤーほど0に近い数字を選択することになるわけです。

下のグラフは、架空データ(n=500)の分布になります。このデータからは33や0、また22付近の数字や、15付近の数字の選択頻度が高いことがわかります。
f:id:hiroecon:20180610092728p:plain



この節では、全ての被験者はLevel-kモデルに基づき数字を選択しているものと仮定し、各レベルの被験者の割合を推定します。推定に際し、次のように仮定します。

  • Level-0の推論を行うプレイヤーは0から100の数字を等確率で選ぶ
  • Level-j(j>0)の推論を行うプレイヤーについては、選択する数字 yは各レベルの最適な推論の数字 y_{j}^{\ast}に平均0、分散 \sigma^{2}正規分布に従う誤差項を加えた数字とする

 y|\text{level}\;j=y_{j}^{\ast}+\epsilon,\;\;\;\;\epsilon\sim N(0,\sigma^{2})

なお、 y_{1}^{\ast}=33 y_{2}^{\ast}=22 y_{3}^{\ast}=15 y_{4}^{\ast}=10 y_{5}^{\ast}=0となります。


これらの仮定によって、各レベルの条件付き確率密度関数
 f(y|\text{level 0})=\frac{1}{100},\;\;0\leq y\leq 100
 f(y|\text{level}\;j)=\frac{1}{\sigma}\phi \Bigl(\frac{y-y_{j}^{\ast}}{\sigma}\Bigr),\;\;0\leq y\leq 100,\;\;j=1,2,3,4,5
となります( \phi(\cdot)正規分布確率密度関数)。


架空データの被験者500名は、Level-0からLevel-5までのプレイヤーで構成されているとし、 p_{j}をLevel-jの被験者の割合とします。これにより、対数尤度関数は次のような形になります。
\displaystyle LogL=\Sigma_{i=1}^{500}\ln \Bigl\{\frac{p_{0}}{100}+\Sigma_{j=1}^{5}p_{j}\frac{1}{\sigma}\phi \Bigl(\frac{y_{i}-y_{j}^{\ast}}{\sigma}\Bigr) \Bigr\}


Rプログラムは以下のようになります。

# read STATA dta file
library(foreign)
d<-read.dta("/Users/Hiro/Documents/R/moffatt/data_files/beauty_sim.dta")
y<-d$y
summary(y)
hist(y,
     breaks=seq(0,100,1),
     xlab="Guess",
     main="Distribution of Simulated Guesses")

negLogL<-function(p1,p2,p3,p4,p5,sigma){
  LogL<-sum(log((1-p1-p2-p3-p4-p5)*(1/100)
                +p1*(1/sigma)*dnorm((y-33)/sigma)
                +p2*(1/sigma)*dnorm((y-22)/sigma)
                +p3*(1/sigma)*dnorm((y-15)/sigma)
                +p4*(1/sigma)*dnorm((y-10)/sigma)
                +p5*(1/sigma)*dnorm((y-0)/sigma)))
  negLogL<--LogL
}

library(stats4)
fit<-mle(negLogL, list(p1=0.3,p2=0.4,p3=0.1,p4=0.1,p5=0.05,sigma=2))
summary(fit)

coef(fit) # estimated coefficients
vcov(fit) # variance-covariance matrix

p1_hat<-coef(fit)[1]
p2_hat<-coef(fit)[2]
p3_hat<-coef(fit)[3]
p4_hat<-coef(fit)[4]
p5_hat<-coef(fit)[5]
sigma_hat<-coef(fit)[6]

p0_hat<-1-p1_hat-p2_hat-p3_hat-p4_hat-p5_hat
p0_hat

# delta method
library(msm)
p0_std.error<-deltamethod(~1-x1-x2-x3-x4-x5,coef(fit),vcov(fit)) # The variables must be labelled x1, x2,...
p0_std.error

# posterior level probabilities
posterior<-function(x) {
  denom<-p0_hat*(1/100)+p1_hat*(1/sigma_hat)*dnorm((x-33)/sigma_hat)+p2_hat*(1/sigma_hat)*dnorm((x-22)/sigma_hat)+p3_hat*(1/sigma_hat)*dnorm((x-15)/sigma_hat)+p4_hat*(1/sigma_hat)*dnorm((x-10)/sigma_hat)+p5_hat*(1/sigma_hat)*dnorm((x-0)/sigma_hat)
  post_0<-p0_hat*(1/100)/denom
  post_1<-p1_hat*(1/sigma_hat)*dnorm((x-33)/sigma_hat)/denom
  post_2<-p2_hat*(1/sigma_hat)*dnorm((x-22)/sigma_hat)/denom
  post_3<-p3_hat*(1/sigma_hat)*dnorm((x-15)/sigma_hat)/denom
  post_4<-p4_hat*(1/sigma_hat)*dnorm((x-10)/sigma_hat)/denom
  post_5<-p5_hat*(1/sigma_hat)*dnorm((x-0)/sigma_hat)/denom
  return(list(post_0,post_1,post_2,post_3,post_4,post_5))
}

guess<-c(1:100)
posterior_type0<-posterior(guess)[[1]]
posterior_type1<-posterior(guess)[[2]]
posterior_type2<-posterior(guess)[[3]]
posterior_type3<-posterior(guess)[[4]]
posterior_type4<-posterior(guess)[[5]]
posterior_type5<-posterior(guess)[[6]]

plot(guess,posterior_type0,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=1)
par(new=T)
plot(guess,posterior_type1,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=2)
par(new=T)
plot(guess,posterior_type2,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=3)
par(new=T)
plot(guess,posterior_type3,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=4)
par(new=T)
plot(guess,posterior_type4,type='l',ann=F,xlim=c(0,100),ylim=c(0,1),lty=5)
par(new=T)
plot(guess,posterior_type5,type='l',xlab="Guess",ylab="Posterior Probability",xlim=c(0,100),ylim=c(0,1),lty=6)
legend_names<-c("posterior type 0","posterior type 1","posterior type 2","posterior type 3","posterior type 4","posterior type 5")
legend("topright",legend=legend_names,lty=c(1:6))


推定結果はこちらになります。

Coefficients:
        Estimate Std. Error
p1    0.39762922 0.02388147
p2    0.11127580 0.01644880
p3    0.08911294 0.01594946
p4    0.04588258 0.01357794
p5    0.04982359 0.01180616
sigma 1.96095481 0.10848212

-2 log L: 3981.022 

なお、Level-0の割合の推定値 \;\hat{p}_{0}は0.3062759、デルタ法でもとめた標準誤差は0.02925361でした。Experimetrics本の411ページの結果とほぼ同じ推定結果が得られました。


最後に、推定値とベイズの定理を用いて、各レベルの事後分布をもとめてみます。例えば、数字 yを観測した条件のもとでの各レベルの事後分布は次のように求めることができます。
\displaystyle f(\text{level}\;0|y)=\frac{  \frac{\hat{p}_{0}}{100}   }{   \frac{\hat{p}_{0}}{100} + \Sigma_{k=1}^{5}\;\frac{1}{\hat{\sigma}}\phi \Bigl(\frac{y_{i}-y_{j}^{\ast}}{\hat{\sigma}}\Bigr)\hat{p}_{k}   }

\displaystyle f(\text{level}\;j|y)=\frac{   \frac{1}{\hat{\sigma}}\phi \Bigl(\frac{y_{i}-y_{j}^{\ast}}{\hat{\sigma}}\Bigr)\hat{p}_{j}   }{   \frac{\hat{p}_{0}}{100} + \Sigma_{k=1}^{5}\;\frac{1}{\hat{\sigma}}\phi \Bigl(\frac{y_{i}-y_{j}^{\ast}}{\hat{\sigma}}\Bigr)\hat{p}_{k}   }\;\;j=1,2,3,4,5


各レベルの事後分布のグラフです。Experimentrics本の412ページのFigure 17.2を再現できました。もし yが40以上ならその被験者はほぼ100%の確率でLevel-0のプレイヤーであることを意味します。他のレベルの事後分布は一つのピークしか無く、レベルが上がるほどピークの yの位置が0に近くなっていくのが特徴です。
f:id:hiroecon:20180610114224p:plain

Experimetrics: 16.4 Estimation of the QRE model

今回はExperimetrics本の16.3.2で求めたpursue-evade gameのQRE(Quantal Response Equilibrium)を使って、架空の被験者データから被験者の合理性を推定します。


推定に際し、次のように仮定します。

  • pursuer役の被験者について、各戦略の選択確率はどの被験者でも同じとする。また各戦略の選択確率は全てのラウンドを通して不変。
  • evaderの役割の被験者についても、各戦略の選択確率はどの被験者でも同じとする。また各戦略の選択確率は全てのラウンドを通して不変。

これにより、合理性を表すエラーパラメーターの推定には、それぞれの役割の各戦略の選択頻度の情報のみが必要となり、対数尤度関数は次のような形になります。
\[
LogL(\mu)=n_{PL}\ln p_{PL}(\mu)+n_{PR}\ln p_{PR}(\mu)+n_{EL}\ln p_{EL}(\mu)+n_{ER}\ln p_{ER}(\mu)
\]


Experimetrics本では、被験者は1ペア(1 pursuer、1 evader)のみで1000ラウンドプレイするとし、各プレイヤーの各戦略の選択頻度は次の値を仮定します。
\[
n_{PL}=500,\;\;\;n_{PR}=500,\;\;\;n_{EL}=667,\;\;\;n_{ER}=333
\]


QREはExperimetrics本の16.3.2で計算したものを使用してます。
hiroecon.hatenablog.com


Rプログラムは以下のようになります。

n1_0<-500 # number of times that pursuer chose L
n1_1<-500 # number of times that pursuer chose R
n2_0<-667 # number of times that evader chose L
n2_1<-333 # number of times that evader chose R

# load mu_p_q (computed in 16.3.2)
mu_p_q<-read.csv("/Users/Hiro/Documents/R/moffatt/mu_p_q.csv")

# log likelihood function
LL<-numeric(nrow(mu_p_q))
for (i in 1:nrow(mu_p_q)) {
  LL[i]<-n1_0*log(mu_p_q$p[i])+n1_1*log(1-mu_p_q$p[i])+n2_0*log(mu_p_q$q[i])+n2_1*log(1-mu_p_q$q[i])
}

LL[which.max(LL)]
mu_p_q$mu[which.max(LL)]
mu_p_q$p[which.max(LL)]
mu_p_q$q[which.max(LL)]

plot(mu_p_q$mu,LL)
abline(v=mu_p_q$mu[which.max(LL)])
text(1.7,-1385,expression(mu==0.72))


横軸は合理性の程度を表すエラーパラメーターの値、縦軸は対数尤度関数の値です。対数尤度関数の値が最大になるのは、エラーパラメーターの値が0.72のときになります。Experimetrics本の結果と同じです。
f:id:hiroecon:20180528085601p:plain

エラーパラメーターの値が0.72のときのQREにおける各戦略の選択確率は次の通り。
\begin{align*}
p_{PL}(\mu)&=0.5001529\\
p_{PR}(\mu)&=0.4998471\\
p_{EL}(\mu)&=0.6668134\\
p_{ER}(\mu)&=0.3331866
\end{align*}

Experimetrics: 16.3.2 Computing the probabilities in the QRE model

今回はExperimetrics本の16.3.2で紹介されているQRE(Quantal Response Equilibrium)の数値計算をRでやってみました。


Rosenthal et al. (2003)に登場するpursue-evadeゲームのQREを求めます。

link.springer.com


ゲームの利得表は次のとおりです。各プレイヤー、0と1の2つの戦略を有します。

Player 2
0 1
Player 1 0 1,-1 0,0
1 0,0 2,-2

元のゲームでは各プレイヤーの戦略はLとRです。ここではLを0、Rを1としています。そして、Player 1(pursuer)が戦略0を選ぶ確率をp、Player 2(evader)が戦略0を選ぶ確率をqと置いています。このとき、Player 1の各戦略の期待利得はそれぞれ
\begin{align*}
ep1_0(q)&=1\times q+0\times (1-q)=q\\
ep1_1(q)&=0\times q+2\times (1-q)=2(1-q)
\end{align*}
となります。同様に、Player 2の各戦略の期待利得はそれぞれ
\begin{align*}
ep2_0(p)&=-1\times p+0\times (1-p)=-p\\
ep2_1(p)&=0\times p+(-2)\times (1-p)=-2(1-p)
\end{align*}
となります。求めたいQREは、
\begin{align*}
p&=\frac{e^{\frac{ep1_0(q)}{\mu}}}{e^{\frac{ep1_0(q)}{\mu}}+e^{\frac{ep1_1(q)}{\mu}}}\\
q&=\frac{e^{\frac{ep2_0(p)}{\mu}}}{e^{\frac{ep2_0(p)}{\mu}}+e^{\frac{ep2_1(p)}{\mu}}}
\end{align*}
という非線形連立方程式を満たすpとqになります。ギリシャ文字のミューはプレイヤーの合理性の程度を表すエラーパラメーターです。0に近づくほど完全合理的になり、無限に近づくほど完全非合理的(ランダムに各戦略を選択)になります。


以下がRのプログラムになります。

# player 1's payoffs
u1_00<-1
u1_01<-0
u1_10<-0
u1_11<-2

# player 2's payoffs
u2_00<--1
u2_10<-0
u2_01<-0
u2_11<--2

fn<-function(x,mu){
  p<-x[1]
  q<-x[2]
  ep1_0<-u1_00*q+u1_01*(1-q)
  ep1_1<-u1_10*q+u1_11*(1-q)
  ep2_0<-u2_00*p+u2_10*(1-p)
  ep2_1<-u2_01*p+u2_11*(1-p)
  belief1<-exp(ep1_0/mu)/(exp(ep1_0/mu)+exp(ep1_1/mu))
  belief2<-exp(ep2_0/mu)/(exp(ep2_0/mu)+exp(ep2_1/mu))
  c(p-belief1,q-belief2)
}

library(nleqslv)
mu_vec<-seq(0.01,10,0.01)
niter<-length(mu_vec)
p_vec<-numeric(niter)
q_vec<-numeric(niter)
for (i in 1:niter) {
  fn_mu<-function(x){
    fn(x,mu_vec[i])
  }
  solutions<-nleqslv(c(0.5, 0.5), fn_mu)
  p_vec[i]<-solutions$x[1]
  q_vec[i]<-solutions$x[2]
}
plot(mu_vec,p_vec,ylim=c(0,1),ylab="probabilities",pch=16,col='red') 
par(new=T)                                          
plot(mu_vec,q_vec,ylim=c(0,1),ylab="",pch=16,col='blue')
abline(h=0.5)
legend("topright", legend = c('p','q'), col = c('red','blue'), pch = c(16,16))

mu_p_q<-data.frame(mu=mu_vec,p=p_vec,q=q_vec)
write.csv(mu_p_q,
          "/Users/Hiro/Documents/R/moffatt/mu_p_q.csv",
          row.names=FALSE)

ここで求めるQREはExperimentrics本16.4で被験者データから合理性の程度を表すエラーパラメーターの推定に用いるため、プログラムの最後でcsvファイルとして書き出しています。
hiroecon.hatenablog.com

プログラムではエラーパラメーターの値を0.01から10まで0.01刻みで変化させ、各ミューの値に対し上記の非線形連立方程式の解であるpとqを求めます。非線形連立方程式を解く必要があるので、プログラムの中でnleqslvというパッケージを使用しています。パッケージの詳細についてはこちらを御覧ください。

CRAN - Package nleqslv


求めたQREを図にしてみました。横軸はエラーパラメータのmuの値(0.01から10まで0.01刻み)、縦軸は確率を表しています。Experimentrics本p.393の図16.2と同じ図が描けました。ミューの値が0に近づくと、pとqの値がそれぞれ混合戦略ナッシュ均衡に収束していきます。一方ミューの値が大きくなるについて、pもqも0.5(つまり、各戦略を等確率で選ぶ)に近づいているのがわかると思います。

f:id:hiroecon:20180526163055p:plain

Experimetrics: 9.3.1 Simulating data from a linear model

今回はExperimetrics本の9.3.1を再現してみます。先に言っておきます。この9.3.1は誤植が多いです。


9.3.1では、まず以下のような線形モデルからシミュレーションデータを生成することになってます(p.214の9.1を参照してください)。
\begin{align*}
y_{i} &= 2.0 + 0.5 x_{i} +u_{i}\;\;\;\;i=1,2,\dots,100\\
x_{i}&\sim U(0,1)\\
u_{i}&\sim N(0,1)
\end{align*}

しかし、実際には
\begin{align*}
y_{i} &= 2.0 + x_{i} +0.5 u_{i}\;\;\;\;i=1,2,\dots,100\\
x_{i}&\sim U(0,1)\\
u_{i}&\sim N(0,1)
\end{align*}
からデータが生成されたようです。

とりあえず、9.1のモデルから生成されたシミュレーションデータを使って推定を行います。使用したRコードは次の通りです。

set.seed(123)
iter<-1000
intercept<-numeric(iter)
slope<-numeric(iter)
sigma<-numeric(iter)

for (i in 1:iter) {
  # DGP (data generating process)
  n<-100
  x<-runif(n,0,1)
  u<-rnorm(n,0,1)
  y<-2+0.5*x+u
  
  fit<-lm(y~x)
  intercept[i]<-coef(fit)[1]
  slope[i]<-coef(fit)[2]
  
  res<-residuals(fit)
  sigma[i]<-sqrt(sum(res^2)/(n-2)) # or sigma[i]<-sigma(fit)
}
library(psych)
describe(intercept)
describe(slope)
describe(sigma)

シミュレーションデータをもとにOLSで推定した切片をinterceptに、傾きをslopeに、誤差項の標準誤差をsigmaに保存します。この作業を1000回繰り返すので、結果的に1000セット分の推定値を得ることになります。

結果はこちら。まずはinterceptから。1000回分の推定値の平均が2.01。データ生成に使用したinterceptが2.0ですからよく推定できてます。

   vars    n mean   sd median trimmed  mad  min max range  skew kurtosis   se
X1    1 1000 2.01 0.19   2.01    2.01 0.19 1.38 2.7  1.32 -0.02     0.14 0.01

次にslope。推定値の平均は0.49。データ生成時のslopeは0.5ですから、こちらも推定できてます。

   vars    n mean   sd median trimmed  mad   min  max range skew kurtosis   se
X1    1 1000 0.49 0.34    0.5    0.49 0.33 -0.67 1.59  2.26 0.02     0.13 0.01

最後にsigma

   vars    n mean   sd median trimmed  mad  min  max range skew kurtosis se
X1    1 1000    1 0.07      1       1 0.07 0.79 1.27  0.49 0.23     0.31  0

1000回分の推定値の分布はこんな感じです。
f:id:hiroecon:20180524083739p:plain


肝心の本の結果の再現ですが、プログラム中のDGPを次のように変更すると、結果を再現できます。

  x<-runif(n,0,1)
  u<-rnorm(n,0,1)
  y<-2+x+0.5*u

Experimetrics by Peter G. Moffatt

今回手に入れた本がこちら。英国・イーストアングリア大学のPeter Moffatt教授の「Experimetrics」です。

Experimetrics: Econometrics for Experimental Economics

Experimetrics: Econometrics for Experimental Economics

サブタイトルにあるように、Experimental Economics(実験経済学)向けのEconometrics(計量経済学)本になります。実験経済学の手法を用いた実験を行なう人にとっては必読書だと思います。

以前の職場では毎年夏に大学院生向けのサマースクールが開講されたのですが、その時にMoffatt教授が今回の本と同じタイトルの講義を1週間の集中講義として開講して下さいました。残念ながら私の能力が低く、そのときはその内容の有用性をあまり理解できておりませんでした。


目次はこちら
1. Introduction and Overview
2. Statistical Aspects of Experimental Design in Experimental Econometrics
3. Treatment Testing
4. Theory Testing, Regression and Dependence
5. Modelling of Decision Times using Regression Analysis
6. Dealing with Discreteness in Experimental Data
7. Ordinal Data in Experimetrics
8. Dealing with Heterogeneity: Finite Mixture Models
9. Simulating Experimental Data, and the Monte-Carlo Method
10. Introduction to the Method of Maximum Simulated Likelihood
11. Dealing with Zeros: Hurdle Models
12. Choice under Risk: Theoretical Issues
13. Choice under Risk: Econometric Modelling
14. Optimal Design in Binary Choice Experiments
15. Social Preference Models
16. Repeated Games and Quantal Response Models
17. Depth of Reasoning Models
18. Learning Models
19. Summary and Conclusion
Appendix A: List of Data Files and Other Files
Appendix B: List of STATA Commands
Appendix C: Choice Problems used in Chapters 5 and 13
References.

サポートページから、使用されるSTATAのソースコードとデータがダウンロード可能です。
www.macmillanihe.com


需要があるかはわかりませんが、自分の勉強も兼ねて、Rを使った再現を今後少しずつやっていこうと思っています。

Rユーザーの方へ:サポートページで提供されているデータはSTATA用のデータですが、Rのforeignパッケージのread.dta()を使ってインポートできます。ちょうどこんな感じです。

# read STATA dta file
library(foreign)
d<-read.dta("/Users/Hiro/Documents/R/moffatt/data files/acquire_sim.dta")

ディレクトリ名は適宜変更してください。