10 Application: Strategy Frequency Estimation

One of the econometric models that most sparked my interest as a grad student was the Strategy Frequency Estimation Method (SFEM), which first showed up in Dal Bó and Fréchette (2011) and Fudenberg, Rand, and Dreber (2012). In these applications, and a lot that followed, the SFEM is applied to data from infinitely repeated Prisoner’s Dilemma experiments. In these games, strategies can be quite complicated beasts, because they must prescribe an action to take after every possible history of play, and histories in an infinitely repeated game can be, well … infinite. We therefore cannot observe strategies (or the pure-action realizations of mixed strategies) like we would in a single-shot game.42 However we cannot just ignore strategies in these games, because a lot of the theory motivating why cooperation will or will not occur has to do with strategies that support cooperation. The idea behind the SFEM is to focus on a menu of strategies that we have decided could be empirically relevant. We then estimate the fraction of participants whose behavior in our experiment is best described by each of these strategies. Hence, the SFEM is a mixture model.

A strategy \(s\) in an indefinitely repeated game can be defined as a function that maps the history of play \(h_{i,t-1}\) into a probability distribution over the actions in player \(i\)’s choice set in the next period of the game, \(y_{i,t}\mid h_{i,t-1},s\), where \(y_{i,t}\) is the action taken by participant \(i\) in period \(t\). While these strategies can in principle be mixed, unless they are fully mixed they will not assign positive probability to all actions, and so we will run into the zero-likelihood problem. The SFEM’s solution to this is to assume that players tremble. That is, participants follow the action specified by the strategy with probability \(1-\epsilon\), and choose another action with probability \(\epsilon\). In the most common application of the SFEM, which is indefinitely repeated Prisoner’s dilemmas, players’ choice sets have just two elements (cooperate and defect), so it makes sense that \(\epsilon\in(0,0.5)\). This ensures that players are more likely to choose the action specified by their strategy than not. For this chapter, I will focus on the SFEM when the choice set is binary43 I will further focus only on strategies that specify pure actions.

The strategy \(s\), tremble probability \(\epsilon\), history \(\{h_{i,t-1}\}\), and data \(\{y_{i,t}\}\) is sufficient information to write down the probability of an individual decision:

\[ p(y_{i,t}\mid h_{i,t-1},s,\epsilon)=(1-\epsilon)^{I(y_{i,t}=s(h_{i,t-1}))}\epsilon^{I(y_{i,t}\neq s(h_{i,t-1}))} \]

where \(I(\cdot)\) is the indicator function. That is, when \(y_{i,t}=s(h_{i,t-1})\) the participant has followed the prescribed action of the strategy, which happens with probability \(1-\epsilon\), and when \(y_{i,t}\neq s(h_{i,t-1})\) the participant has not followed the prescribed action of the strategy, which happens with probability \(\epsilon\).

As the goal of the SFEM is to estimate the fraction of participants who use each strategy, we combine these decision-level likelihoods and mixing probabilities \(\{\rho_s\}_{s=1}^S\) into a grand likelihood using a dichotomous mixture model specification:

\[ \begin{aligned} \log p(y\mid \rho,\epsilon) &=\sum_{i=1}^N\log\left(\sum_{s=1}^S\exp\left(\log \rho_s+\sum_{t=1}^{T_i}\log p(y_{i,t}\mid h_{i,t-1},s,\epsilon)\right)\right) \end{aligned} \]

This is one of these situations where we can take advantage of Stan’s log_sum_exp function, so I have added the logged mixing probabilities to the exponent, rather than multiplying them in levels. Also, it is important to keep track of the \(i\) subscript on \(T_i\), the total number of actions taken by participant \(i\). Indefinitely repeated games is one of the situations where we do not have a balanced panel, as the length of a game must be random.

10.1 Simplifying the individual likelihood functions

For pure strategies, there is a lot we can do to speed up estimation. Specifically, we can simplify the summation over \(t\) in the above grand likelihood by pre-summarizing our data: we don’t actually need to know for which specific decisions a participant followed a strategy or trembled, we just need to know the total number of times they followed a strategy, and the total number of times they trembled. This is because we have assumed that actions are independent conditional on \(h_{i,t-1}\), \(s\), and \(\epsilon\). To see this, note that:

\[ \begin{aligned} \sum_{t=1}^{T_i}\log p(y_{i,t}\mid h_{i,t-1},s,\epsilon) &=\sum_{t=1}^{T_i}\left[I(y_{i,t}=s(h_{i,t-1}))\log(1-\epsilon)+I(y_{i,t}\neq s(h_{i,t-1}))\log (\epsilon)\right]\\ &=\log(1-\epsilon)\sum_{t=1}^{T_i}I(y_{i,t}=s(h_{i,t-1}))+\log(\epsilon)\sum_{t=1}^{T_i}I(y_{i,t}\neq s(h_{i,t-1}))\\ &=\log(1-\epsilon)\mathrm{follow}_{i,s}+\log(\epsilon)(T_i-\mathrm{follow}_{i,s}) \end{aligned} \]

where \(\mathrm{follow}_{i,s}\) is the number of times participant \(i\) followed the prescribed action in strategy \(s\). Importantly, \(\mathrm{follow}_{i,s}\) can be pre-calculated, so we can speed up our estimation by doing this.

10.2 Example experiment: Dal Bó and Fréchette (2011)

In their Table 7, Dal Bó and Fréchette (2011) estimate the fraction of \(S=6\) strategies being used in their indefinitely repeated Prisoner’s dilemma experiment.44 These strategies are:

  1. Always Defect (AD)
  2. Always Cooperate (AC)
  3. Grim Trigger (G), defect for ever once one person has defected
  4. Tit for Tat (TFT)
  5. Win Stay Lose Shit (WSLS)
  6. Tit for two tats (T2)

From their footnote 19:

WSLS is a strategty hat starts cooperating and then condition behavior only on the outcome of the previous round. If either both cooperate or neither cooperate, then WSLS cooperates; otherwise it defects. T2 starts cooperating, and a defection by the other triggers two rounds of defection, after which the strategy goes back to cooperation. These two strategies are cooperative st rategies with punishments of limited length.

Each strategy specifies an action (for an infinitely repeated Prisoner’s Dilemma, either “Cooperate” or “Defect”) to be played after every possible combination of strategies. The trouble, with a few notable exceptions, is that we do not observe strategies, we only observe the actions that result from these strategies, and so we have data like this:

DBFData<-readRDS('Data/SFEM/DBFData.rds')
DBFData |> head() |> kbl() |> kable_classic(full_width=F)
id follow1 follow2 follow3 follow4 follow5 follow6 n Treatment
1 21 4 13 13 8 9 25 1
2 25 0 17 16 1 13 25 1
3 25 0 17 17 0 13 25 1
4 22 3 14 15 3 12 25 1
5 25 0 17 15 2 13 25 1
6 13 12 5 11 11 9 25 1

Each row here is a participant in the experiment. Follow1 is a count of the number of times this participant played the action consistent with AD, likewise for follow2, and so on. \(n\) is the number of actions this participant played (some data from early rounds were excluded).

10.2.1 The SFEM with homogeneous trembles

Here is my implementation of the model in Stan. I chose not to use the log_sum_exp function here because it would have meant that I could not vectorize the calculation of the target += log(like_i*mix) line. In order to make the estimates directly comparable to Table 7 of Dal Bó and Fréchette (2011), I estimate one model for each of their six treatments.45

// SFEM.stan
functions {
   
              
}
data {
  // Tell Stan how to read in the data
  int N; // Number of subjects
  int S; // Number of strategies
  matrix[N,S] COUNTS; // Number of decisions by each player, here I repeat the vector S times to make the stuff below a bit simpler
  matrix[N,S] FOLLOWS; // counts of following each strategy
  vector[S] priorMix; // Dirichlet prior for mixing probabilities
}
transformed data {
    
}
parameters {
  real<lower=0,upper=0.5> trmb; // tremble probability
  simplex[S] mix; // strategy mixing probabilities
  

}
transformed parameters {
  
}
model {
    matrix[N,S] like_i; 
    
    mix ~ dirichlet(priorMix);

  // Individual likelihoods
    like_i = exp(FOLLOWS*log(1-trmb)+(COUNTS-FOLLOWS)*log(trmb));
    
    // Adding the grand likelihood to the target
    target += log(like_i*mix);
    
    


  }
generated quantities {

}

The estimates from these models are shown in Table @(tab:SFEMSimpleTable), which is arranged in almost the same way as Table 7 in Dal Bó and Fréchette (2011).

TAB<-readRDS("Code/SFEM/SFEMHomogeneousTremblesEstimates.rds") 
rownames(TAB)<-c("AD","AC","G","TFT","WSLS","T2","tremble")

TAB |> 
  kbl(caption = "SFEM estimates from the six treatments of @DBF2011. Posterior means with standard deviations in parentheses. Each colum is one of the six treatments. $R$ is the payoff when both payers cooperate, and $\\delta$ is the continuation probability.") |> 
  kable_classic(full_width=F) |>
  add_header_above(c(" "=1,"R = 32"=1,"R = 40"=1,"R = 48"=1,"R = 32"=1,"R = 40"=1,"R = 48"=1)) |>
  add_header_above(c(" "=1,"$\\delta=\\frac{1}{2}$"=3,"$\\delta=\\frac{3}{4}$"=3))
\(\delta=\frac{1}{2}\)
\(\delta=\frac{3}{4}\)
R = 32
R = 40
R = 48
R = 32
R = 40
R = 48
Table 10.1: SFEM estimates from the six treatments of Dal Bó and Fréchette (2011). Posterior means with standard deviations in parentheses. Each colum is one of the six treatments. \(R\) is the payoff when both payers cooperate, and \(\delta\) is the continuation probability.
AD 0.83 (0.053) 0.716 (0.061) 0.49 (0.07) 0.59 (0.071) 0.118 (0.048) 0.02 (0.021)
AC 0.02 (0.02) 0.083 (0.04) 0.084 (0.04) 0.021 (0.021) 0.27 (0.094) 0.112 (0.072)
G 0.02 (0.019) 0.051 (0.032) 0.035 (0.033) 0.034 (0.029) 0.225 (0.096) 0.145 (0.1)
TFT 0.09 (0.042) 0.107 (0.045) 0.332 (0.071) 0.314 (0.069) 0.297 (0.11) 0.455 (0.117)
WSLS 0.02 (0.019) 0.021 (0.02) 0.037 (0.027) 0.02 (0.019) 0.045 (0.043) 0.057 (0.053)
T2 0.02 (0.019) 0.022 (0.022) 0.022 (0.022) 0.02 (0.02) 0.045 (0.041) 0.211 (0.112)
tremble 0.06 (0.008) 0.137 (0.01) 0.089 (0.008) 0.097 (0.007) 0.092 (0.01) 0.03 (0.005)

Comparing my estimates in Table 10.1 to Table 7 of Dal Bó and Fréchette (2011), there are many similarities, but my mixing probabilities are all closer to \(\frac{1}{6}\) (i.e. large probabilities are attenuated, small probabilities are inflated). This is because \(\frac{1}{6}\) is the prior mean, so we should expect this. Some common themes between these tables are:

  • Mostly AD in the \(R=32\) or \(R=40\), \(\delta=\frac12\) treatments
  • Similar fractions of AD and TFT in the \(R=48\) \(\delta=\frac12\) and \(R=32\) \(\delta=\frac34\) treatments
  • Similar fractions of TFT in the \(\delta=\frac34\) \(R=40\) and \(R=48\) treatments, and
  • Similar fraction of AC in the \(\delta=\frac34\) \(R=40\) treatment.

10.2.2 Adding heterogeneous trembles and integrating the likelihood

An extension I did with the SFEM in Bland (2020) was to allow for participant-specific heterogeneous tremble probabilities \(\epsilon_i\). This can be done a few different ways. Probably the most obvious way is to assume that \(2\epsilon_i\) probit-normally distributed,46 and use data augmentation accordingly. I use this distribution in Bland (2020), but implement a maximum likelihood estimator that uses Monte Carlo integration to integrate out the \(\epsilon_i\)s.47 Since the data augmentation approach is going to look like any other 1-parameter hierarchical extension that I might show you in this book, I thought I would show you a different solution, where we can find an analytical solution to integrating the likelihood. This is similar to the process of finding conjugate priors.

To begin with, I will write out the strategy-specific likelihood function for one participant’s decisions if we assume a general probability density function for the distribution of \(\epsilon\), call it \(p(\epsilon\mid \theta)\), where \(\theta\) are some parameters describing the distribution.

\[ \begin{aligned} p(y_i\mid h_{i,t-1},s,\theta)&=\int_0^{0.5}\prod_{t=1}^{T_i}p(y_{i,t}\mid h_{i,t-1},s,\epsilon_i)p(\epsilon\mid \theta)\mathrm d\epsilon\\ &=\int_0^{0.5}(1-\epsilon)^{F_{i,s}}\epsilon^{T_i-F_{i,s}}p(\epsilon\mid\theta)\mathrm d\epsilon \end{aligned} \]

where \(F_{i,s}\) is short for the variable \(\mathrm{follow}_{i,s}\).

Looking at the likelihood part, it almost looks like an incomplete beta function, which is:

\[ B(x;a,b)=\int_0^xt^{a-1}(1-t)^{b-1}\mathrm dt \]

so if we choose:

\[ p(\epsilon\mid \theta)\propto \epsilon^{\theta_1-1}(1-\epsilon)^{\theta_2-1} \]

our likelihood becomes:

\[ \begin{aligned} p(y_i\mid h_{i,t-1},s,\theta)&\propto\int_0^{0.5}(1-\epsilon)^{F_{i,s}}\epsilon^{T_i-F_{i,s}}\epsilon^{\theta_1-1}(1-\epsilon)^{\theta_2-1}\mathrm d\epsilon\\ &=\int_0^{0.5}\epsilon^{\theta_1+T_i-F_{i,s}-1}(1-\epsilon)^{\theta_2+F_{i,s}-1}\mathrm d\epsilon\\ &=B(0.5;\theta_1+T_i-F_{i,s},\theta_2+F_{i,s}) \end{aligned} \]

So we can evaluate the likelihood, up to a constant of proportionality, using the incomplete Beta function! Now we need to work out the constant of proportionality. Fortunately for us, this is just the constant to make sure that \(p(\epsilon\mid\theta)\) integrates to one, so it must be that:

\[ \begin{aligned} 1&=c\int_0^{0.5}\epsilon^{\theta_1-1}(1-\epsilon)^{\theta_2-1}\mathrm d\epsilon\\ &=cB(0.5;\theta_1,\theta_2)\\ c&=\frac{1}{B(0.5;\theta_1,\theta_2)}\\ \implies p(\epsilon\mid\theta)&=\frac{\epsilon^{\theta_1-1}(1-\epsilon)^{\theta_2-1}}{B(0.5;\theta_1,\theta_2)} \end{aligned} \]

And so we can write the strategy-specific log-likelihood of a participant’s choices as:

\[ \begin{aligned} \log p(y_i\mid h_{i,t-1},s,\theta)&=\log B(0.5;\theta_1+T_i-F_{i,s},\theta_2+F_{i,s})-\log B(0.5; \theta_1,\theta_2) \end{aligned} \]

This means that we have modeled the tremble probabilities \(\epsilon_i\) as coming from a truncated beta distribution.

Unfortunately, Stan’s inc_beta function returns the regularized incomplete beta function, which is not what we have above. The regularized version divides the incomplete beta function by the (complete) beta function. That is, we can evaluate:

\[ \begin{aligned} \text{(complete) beta: } B(a,b)&=\int_0^1t^{a-1}(1-t)^{b-1}\mathrm dt\\ \text{regularized incomplete beta: } I_x(a,b)&=\frac{\int_0^xt^{a-1}(1-t)^{b-1}\mathrm dt}{B(a,b)} \end{aligned} \]

so we will need to add in the beta function to this expression. In summary, we have:

\[ \begin{aligned} \log p(y_i\mid h_{i,t-1},s,\theta)&=\log I_{0.5}(\theta_1+T_i-F_{i,s},\theta_2+F_{i,s})+\log B(\theta_1+T_i-F_{i,s},\theta_2+F_{i,s})\\ &\quad -\log I_{0.5}(\theta_1,\theta_2)-\log B(\theta_1,\theta_2) \end{aligned} \]

In Stan inc_beta is the regularized incomplete beta function, and lbeta is the natural logarithm of the (complete) beta function.

Now that we have two more parameters, \(\theta_1\) and \(\theta_2\), we need to assign priors for these. When using the (un-truncated) beta distribution, I like to re-parameterize this into a mean parameter and a strength parameter like this:

\[ \begin{aligned} q&=\frac{\theta_1}{\theta_1+\theta_2}\in(0,1)\\ k&=\theta_1+\theta_2>0 \end{aligned} \]

for the un-truncated beta distribution, we can therefore think about it as having a mean \(q\) and concentration \(k\). The interpretation is less clear for this truncated case, but I will go ahead with it and use this re-parametrized version to select priors. Since \(q\) is the mean, and can take on any probability between 0 and 1,48 I choose a uniform prior. For \(k\), I use a fairly spread-out:

\[ k\sim \mathrm{Cauchy}^+(0,1) \]

Here is the code I came up with to implement this model in Stan. I was somewhat disappointed that Stan didn’t have a vectorized version of the inc_beta function, so I had to drop my vectorized implementation of calculating like_i from the homogeneous model. All of this is to say, I am not sure whether this implementation where we analytically solve for the integrated likelihood is actually faster than using data augmentation. This is because with data augmentation I would be able to vectorize the operations.

// SFEMBeta.stan
functions {
   
              
}
data {
  // Tell Stan how to read in the data
  int N; // Number of subjects
  int S; // Number of strategies
  matrix[N,S] COUNTS; // Number of decisions by each player, here I repeat the vector S times to make the stuff below a bit simpler
  matrix[N,S] FOLLOWS; // counts of following each strategy
  vector[S] priorMix; // Dirichlet prior for mixing probabilities
  
  real<lower=0> prior_k;
}
transformed data {
    
}
parameters {
  
  simplex[S] mix; // strategy mixing probabilities
  
  
  real<lower=0,upper=1> q; // mode of tremble distribution
  real<lower=0> k; // concentration of tremble distribution

}
transformed parameters {
  
  vector[2] theta;
  
  theta[1] = k*q;
  theta[2] = k*(1-q);
  
}
model {
    matrix[N,S] like_i; 
    
    mix ~ dirichlet(priorMix);
  
  // Individual likelihoods
  
  /* As far as I can tell, there is not a vectorized version of 
  inc_beta. Hence the double for loop below here.
  
  Another thing that tripped me up was that inc_beta takes arguments
  inc_beta(a,b,x), not inc_beta(x,a,b). Always read the manual, peeps!
  */
  for (ii in 1:N) {
    for (ss in 1:S) {
      like_i[ii,ss] = exp(
                  log(inc_beta(theta[1]+COUNTS[ii,ss]-FOLLOWS[ii,ss],theta[2]+FOLLOWS[ii,ss],0.5))
                  +lbeta(theta[1]+COUNTS[ii,ss]-FOLLOWS[ii,ss],theta[2]+FOLLOWS[ii,ss])
                  -log(inc_beta(theta[1],theta[2],0.5))
                  -lbeta(theta[1],theta[2])
                );
    }
  }

    // Adding the grand likelihood to the target
    target += log(like_i*mix);
    
    // prior for k
    k~ cauchy(0.0,prior_k);

  }
generated quantities {

}

Similarly to the previous section, I estimate one model for each treatment in Dal Bó and Fréchette (2011). The estimates are summarized in Table 10.2. These estimates are very similar to those in Table 10.1, and therefore we would draw similar conclusions about strategy frequencies with these heterogeneous-trembles models. In this case the heterogeneity does not seem to be changing the results that we are most interested in much.

TAB<-readRDS("Code/SFEM/SFEMHeterogeneousTremblesEstimates.rds") 
rownames(TAB)<-c("AD","AC","G","TFT","WSLS","T2","q","k")

TAB |> 
  kbl(caption = "SFEM estimates from the six treatments of @DBF2011, assuming a heterogeneous distribution of tremble probability $\\epsilon$. Posterior means with standard deviations in parentheses. Each colum is one of the six treatments. $R$ is the payoff when both payers cooperate, and $\\delta$ is the continuation probability.") |> 
  kable_classic(full_width=F) |>
  add_header_above(c(" "=1,"R = 32"=1,"R = 40"=1,"R = 48"=1,"R = 32"=1,"R = 40"=1,"R = 48"=1)) |>
  add_header_above(c(" "=1,"$\\delta=\\frac{1}{2}$"=3,"$\\delta=\\frac{3}{4}$"=3))
\(\delta=\frac{1}{2}\)
\(\delta=\frac{3}{4}\)
R = 32
R = 40
R = 48
R = 32
R = 40
R = 48
Table 10.2: SFEM estimates from the six treatments of Dal Bó and Fréchette (2011), assuming a heterogeneous distribution of tremble probability \(\epsilon\). Posterior means with standard deviations in parentheses. Each colum is one of the six treatments. \(R\) is the payoff when both payers cooperate, and \(\delta\) is the continuation probability.
AD 0.851 (0.051) 0.749 (0.063) 0.481 (0.072) 0.587 (0.071) 0.121 (0.049) 0.02 (0.02)
AC 0.021 (0.02) 0.047 (0.032) 0.086 (0.041) 0.021 (0.02) 0.254 (0.091) 0.11 (0.074)
G 0.021 (0.021) 0.058 (0.034) 0.039 (0.036) 0.047 (0.033) 0.247 (0.094) 0.13 (0.1)
TFT 0.065 (0.037) 0.104 (0.046) 0.337 (0.074) 0.304 (0.068) 0.285 (0.109) 0.471 (0.12)
WSLS 0.021 (0.02) 0.021 (0.02) 0.033 (0.027) 0.02 (0.02) 0.047 (0.044) 0.062 (0.057)
T2 0.021 (0.02) 0.021 (0.021) 0.024 (0.023) 0.021 (0.02) 0.046 (0.044) 0.207 (0.113)
q 0.354 (0.254) 0.547 (0.236) 0.451 (0.245) 0.4 (0.239) 0.389 (0.243) 0.226 (0.234)
k 0.728 (0.814) 0.831 (0.589) 0.676 (0.57) 1.278 (1.146) 0.673 (0.672) 1.461 (1.826)

10.3 R code to do these estimations

library(tidyverse)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

rounding<-3

####################################################
# Loading the data from the files available on the 
# journal website
####################################################


DBFData<-data.frame()
for (ff in 1:6) {
  fname<-paste("Data/SFEM/dfformatlab_strg_",ff,"_special.txt",sep="")
  d<-data.frame(read.delim(fname, header = FALSE, sep = "\t", dec = "."))
  colnames(d)<-c("match","round","treatment","coop","id","ocoop_all","strg1","strg2","strg3","strg4","strg5","strg6","session","id2_all")
  DBFData<-rbind(DBFData,d)
}

for (ff in 1:6) {
  str<-paste("DBFData$follow",ff,"<- as.integer(DBFData$coop == DBFData$strg",ff,")",sep="")
  eval(parse(text=str))
}
AggData<-aggregate(DBFData[,c("follow1","follow2","follow3","follow4","follow5","follow6")], by=list(DBFData$id),FUN=sum)
colnames(AggData)[1]<-"id"
AggData$n<-aggregate(DBFData$id,by=list(DBFData$id),FUN=length)$x
AggData$Treatment<-aggregate(DBFData$treatment,by=list(DBFData$id),FUN=mean)$x
saveRDS(AggData,"Data/SFEM/DBFData.rds")

priorMix<-c(1,1,1,1,1,1)

##############################
# Homogeneous trembles models
##############################

model<-stan_model("Code/SFEM/SFEM.stan")

Fits<-c()

for (tt in 1:6) {
  print(paste("estimating homogeneous trembles model",tt))
  Dt<- AggData %>% filter(Treatment==tt)
  d = list(
    FOLLOWS=Dt[,2:7],
    N=dim(Dt)[1],
    S=6,
    COUNTS=kronecker(Dt[,"n"],matrix(1,1,6)),
    priorMix=priorMix)
  Fit<-sampling(model,data=d,seed=42,control = list(adapt_delta = 0.8))
  S<-summary(Fit,pars=c("mix","trmb"))$summary[,c("mean","sd")] |>
    data.frame() |>
    mutate(msd = paste0(mean |> round(rounding)," (",sd |> round(rounding),")"))
  Fits<-cbind(Fits,S$msd)
}

Fits |> saveRDS("Code/SFEM/SFEMHomogeneousTremblesEstimates.rds")

##############################
# Heterogeneous trembles models
##############################


model<-stan_model("Code/SFEM/SFEMBeta.stan")

Fits<-c()

for (tt in 1:6) {
  print(paste("estimating heterogeneous trembles model",tt))
  Dt<- AggData %>% filter(Treatment==tt)
  d = list(
    FOLLOWS=Dt[,2:7],
    N=dim(Dt)[1],
    S=6,
    COUNTS=kronecker(Dt[,"n"],matrix(1,1,6)),
    priorMix=priorMix,
    prior_k = 1)
  Fit<-sampling(model,data=d,seed=42,control = list(adapt_delta = 0.8))
  S<-summary(Fit,pars=c("mix","q","k"))$summary[,c("mean","sd")] |>
    data.frame() |>
    mutate(msd = paste0(mean |> round(rounding)," (",sd |> round(rounding),")"))
  Fits<-cbind(Fits,S$msd)
}

Fits |> saveRDS("Code/SFEM/SFEMHeterogeneousTremblesEstimates.rds")

References

———. 2020. “Heterogeneous Trembles and Model Selection in the Strategy Frequency Estimation Method.” Journal of the Economic Science Association 6 (2): 113–24.
Dal Bó, Pedro, and Guillaume R Fréchette. 2011. “The Evolution of Cooperation in Infinitely Repeated Games: Experimental Evidence.” American Economic Review 101 (1): 411–29.
———. 2019. “Strategy Choice in the Infinitely Repeated Prisoner’s Dilemma.” American Economic Review 109 (11): 3929–52.
Fudenberg, Drew, David G Rand, and Anna Dreber. 2012. “Slow to Anger and Fast to Forgive: Cooperation in an Uncertain World.” American Economic Review 102 (2): 720–49.
Romero, Julian, and Yaroslav Rosokha. 2018. “Constructing Strategies in the Indefinitely Repeated Prisoner’s Dilemma Game.” European Economic Review 104: 185–219.

  1. Two exceptions to this are Romero and Rosokha (2018) and Dal Bó and Fréchette (2019). ↩︎

  2. For choice sets with more than two elements we need to be a bit more careful with how we specify the trembles. Since I am focusing attention on indefinitely-repeated Prisoner’s dilemma experiments, where the choice set is binary, I will not go into more detail here, but might later in a follow-up chapter.↩︎

  3. You can download their data (and Matlab code to replicate it!) from here↩︎

  4. These estimations run very quickly (a second or two on my laptop each), so it’s probably not worth trying to find a korok while it runs.↩︎

  5. Multiplying by 2 ensures that \(\epsilon_i\in(0,0.5)\).↩︎

  6. This was not for my lack of trying. I originally wrote that paper using Bayesian estimators, but the editor wanted me to re-do everything using MLE. ↩︎

  7. For the truncated model, if \(q<0.5\) then \(q\) is the mode of the distribution.↩︎