Preliminaries

First we load the required libraries. The library used for the empirical game-theoretic analysis is empiricalGameTheory.

require(knitr)
require(fBasics)
require(moments)
require(empiricalGameTheory)

Evolutionary Game-Theory

In an evolutionary game-theoretic model, pairs of agents are chosen at random from an idealised infinite population. By assuming an infinite population and no mutation we can specify the dynamics of the system using a simple ordinary differential equation. The standard approach is to use the replicator dynamics equation (Weibull, 1997) to model how the frequency of each strategy in the larger population changes over time in response to the within-group payoffs:

\(\dot{m}_i = \left[u(e_i,\vec{m}) - u(\vec{m},\vec{m})\right]m_i\)

where \(\vec{m}\) is a mixed-strategy vector, \(u(\vec{m},\vec{m})\) is the mean payoff when all players play \(\vec{m}\), and \(u(e_i,\vec{m})\) is the average payoff to pure strategy \(i\) when all players play \(\vec{m}\), and \(\dot{m}_i\) is the first derivative of \(m_i\) with respect to time. Strategies that gain above-average payoff become more likely to be played, and this equation models a simple co-evolutionary process of adaptation.

Heuristic Payoff Matrices

When computing payoffs in the above equation, the empirical game-theory method uses a payoff matrix which is calibrated by running very many simulations. The payoff matrix is said to be heuristic because several simplifying assumptions are made in the interests of tractability. We can make one important simplification by assuming that the game is symmetric, and therefore that the payoff to a given strategy depends only on the of agents within the group adopting each strategy. Thus for a game with \(j\) strategies, we represent the payoff matrix as a pair \((\vec{p}, \vec{q})\). The vector \(\vec{p} = \left( p_1, \ldots, p_j \right)\) represents the group composition, where \(p_i\) specifies the number of agents who are playing the \(i^{\mathrm{th}}\) strategy. Each entry \(\vec{p} \in P\) is mapped onto an outcome vector \(\vec{q} \in Q\) of the form \(\vec{q} = \left( q_1, \ldots, q_j \right)\) where \(q_i\) specifies the expected payoff to the \(i^{\mathrm{th}}\) strategy.

For a game with \(n\) agents, the number of entries in the payoff matrix is given by $ s = $. For example, for \(n=10\) agents and \(j=5\) strategies, we have a payoff matrix with \(s = 1001\) entries.

For each entry in the payoff matrix we estimate the expected payoff to each strategy by running \(5 \times 10^5\) simulations and taking the mean fitness rounded according to the corresponding standard error.

For example, if we have an entry in the payoff matrix \(\vec{p} = \left( 5, 4, 1, 0, 0 \right)\) then we would run \(10^5\) simulations with \(n=5+4+1+0+0=10\) agents, five of which would be initialised to use the first strategy, four to use the second strategy, one using the third strategy, and zero agents using the remaining strategies.

Paper-Rock-Scissors Example

Initially, rather than using simulations to populate the payoff matrix we will consider a simple example where the payoffs are known apriori without having to resort to simulation.

The game of paper-rock-scissors is a standard normal-form game with 2-players and 3-strategies, labelled rock (\(R\)), paper (\(P\)) and scissors (\(S\)). If we write the number of players adopting a particular strategy as \(N_X\), and the payoff to the player(s) adopting strategy \(X\) as \(u(e_X)\), then the payoff matrix for this game can be specified according to the table below.

data(payoff.matrix.rps)
kable(payoff.matrix.rps, col.names = c('$N_R$', '$N_P$', '$N_S$', '$u(e_R)$', '$u(e_P)$', '$u(e_S)$'))
\(N_R\) \(N_P\) \(N_S\) \(u(e_R)\) \(u(e_P)\) \(u(e_S)\)
0 0 2 0 0 0
0 1 1 0 -1 1
2 0 0 0 0 0
1 1 0 -1 1 0
0 2 0 0 0 0
1 0 1 1 0 -1

The function HeuristicPayoff can be used to estimate the function \(u(e_i, \vec{m})\). For example, suppose we want the payoff to pure strategy \(e_R\) against the strategy \(e_S\), we can express paper as a mixed-strategy where the probability of playing paper is \(1\), i.e. \(\vec{m} = (0, 1, 0)\). In R:

R = 1
P = 2
S = 3
HeuristicPayoff(payoff.matrix.rps, R, c(0, 1, 0))
## [1] -1

We can generalise this to expected payoffs when the second strategy is a vector representing the probability with which each pure strategy is adopted:

So the payoff to rock, when our opponent plays rock with probability \(\frac{1}{5}\), paper with probability \(\frac{1}{5}\) and scissors with probability \(\frac{3}{5}\) can be written as \(u(e_R, (0.2, 0.2, 0.6))\), and in R:

HeuristicPayoff(payoff.matrix.rps, R, c(0.2, 0.2, 0.6))
## [1] 0.4

We can encapsulate this inside an R object belonging to the class HeuristicGame:

game.rps <- HeuristicGameFromPayoffMatrix(payoff.matrix.rps, strategies = c('R', 'P', 'S'))

We can also solve the replicator-dynamics diffential equation using finite-difference methods implementing in R’s deSolve library. The method ComputeTrajectory will take an intial condition in the form of a mixed-strategy representing the initial fraction of the population playing each pure strategy, and then iteratively integrate the equation. So if we start with \(x_0 = (0.2, 0.2, 0.6)\), we can see how the population will evolve as a time-series:

t <- seq(0, 100, by=0.1)
population.frequencies <- as.ts(ComputeTrajectory(game.rps, c(0.2, 0.2, 0.6), 
times = t))
ts.plot(population.frequencies, lty=1:3)

If we do this for a large sample of randomly-chosen initial conditions, then we can plot the phase-space of the corresponding trajectories in a 2-dimensional space:

# Compute trajectories for all specified initial conditions over the specified time interval
if (!var.exists('game.rps.analysed'))
  game.rps.analysed <- Analyse(game.rps, initial.values = initial.values.random, times = t, parallel=T)
result <- plot(game.rps.analysed)

Using simulation to estimate payoffs

Sometimes the payoff matrix is not known in advance, and we may need to use simulation in order to estimate the payoffs to each behavioural rule or strategy. This is called empirical game-theory.

Here we analyse an adaptive-expectations model of a financial market, which incoprorates market micro-structure (Chiarella and Iori, 2002). The payoffs in this example have already been obtained by running the JASA simulator.

In this case we have \(n=12\) agents and three strategies: chartist (\(C\)), fundamentalist (\(F\)) and noise-trader (\(N\)). For a particular profile of strategies, we can run the simulation and observe the resulting payoffs to each strategy.

Here are the payoffs for a single run of the simulation when we have four agents playing each of the three strategies, i.e. \(\vec{m} = ( \frac{1}{3}, \frac{1}{3}, \frac{1}{3} )\).

payoffs <- read.csv('data/payoffs1.csv', sep='\t', header=F)
strategies <- c('C', 'F', 'N')
names(payoffs) <- c('t', strategies)
plot.ts(payoffs[strategies], main='Payoffs to Noise (N), Fundamentalist (F) and Chartist (C)')

In this case, on average the fundamentalist strategy performs marginally better than the noise-trader strategy,

kable(t(sapply(payoffs[strategies], mean)), col.names = strategies)
C F N
0.0431878 0.7528448 -0.3436396

but notice that there is considerable variance and skewness in the outcomes

p <- payoffs[strategies]
Reduce(rbind, Map(function(x) sapply(p, x), c(range, var, skewness, kurtosis)))
##                C         F           N
## [1,]  -8.9056836  0.000000 -230.159946
## [2,]  30.9606588 97.735787  163.919387
## [3,]   0.7735191 23.077706   92.144318
## [4,]   9.6334583  7.800003   -2.007374
## [5,] 197.4800323 78.389310  124.365939
rownames(stats) <- c('min', 'max', 'variance', 'skewness', 'kurtosis')
kable(stats, col.names = strategies)
C F N
min -8.9056836 0.000000 -230.159946
max 30.9606588 97.735787 163.919387
variance 0.7735191 23.077706 92.144318
skewness 9.6333138 7.799886 -2.007344
kurtosis 194.4760827 75.387742 121.363451 
hist(payoffs[,'N'])

hist(payoffs[,'F'])

hist(payoffs[,'C'])

This suggests that risk-sensitivity will be important to the analysis. Also note that the fundamentalist always makes a profit, but noise-traders sometimes make a loss.

Note, however, that the expected payoffs to each strategy will change depending on the profile of strategies in play. This means we need to run simulations for every partition of strategies over agents. Below are the payoffs obtained to each strategy in simulation averaged over 1000 simulation runs for each row when we use a risk-neutral utility function.

columns <- c('$N_C$', '$N_F$', '$N_N$', '$u(e_F)$', '$u(e_C)$', '$u(e_N)$')
game.abm.rn <- HeuristicGame('data/payoff-matrix.csv', strategies, 1, Utility=RiskNeutralUtility)
kable(game.abm.rn$payoffs.normalised, col.names=columns)
\(N_C\) \(N_F\) \(N_N\) \(u(e_F)\) \(u(e_C)\) \(u(e_N)\)
0 0 12 NaN NaN 5.283930
0 1 11 NaN -3.313542 -0.075599
1 0 11 2.684824 NaN 57.338964
0 2 10 NaN 4.392457 -0.005418
1 1 10 0.784595 -0.617543 0.023029
2 0 10 4.213074 NaN 15.512478
0 3 9 NaN 4.418209 -0.037730
1 2 9 0.567423 2.516196 0.039516
2 1 9 0.743009 -1.265723 0.225008
3 0 9 100.000000 NaN -100.000000
0 4 8 NaN 3.398435 -0.022088
1 3 8 0.370492 2.767608 -0.053961
2 2 8 0.376618 2.104009 0.006200
3 1 8 0.453474 -0.672825 -0.016591
4 0 8 2.494006 NaN -15.959641
0 5 7 NaN 2.462643 0.009682
1 4 7 0.238200 2.156917 0.034972
2 3 7 0.238020 2.205446 0.003783
3 2 7 0.300584 1.880427 0.042752
4 1 7 0.274598 0.216531 0.043825
5 0 7 91.798089 NaN -100.000000
0 6 6 NaN 1.742675 -0.028795
1 5 6 0.175443 1.549591 0.118636
2 4 6 0.233682 1.808741 -0.006925
3 3 6 0.183666 1.678503 0.071415
4 2 6 0.179918 1.432854 -0.110305
5 1 6 0.226900 0.290152 0.151204
6 0 6 100.000000 NaN 100.000000
0 7 5 NaN 1.200824 -0.024103
1 6 5 0.170945 1.162645 0.015280
2 5 5 0.121578 1.209381 -0.097011
3 4 5 0.133677 1.330242 0.085419
4 3 5 0.118462 1.297028 0.056280
5 2 5 0.129434 1.299713 0.053112
6 1 5 0.157284 0.497796 0.066019
7 0 5 0.895871 NaN -14.841985
0 8 4 NaN 0.838704 0.010219
1 7 4 0.094206 0.736950 -0.010892
2 6 4 0.092014 0.815534 -0.045597
3 5 4 0.091360 0.882920 -0.001054
4 4 4 0.089292 0.957566 0.000993
5 3 4 0.077106 0.994239 0.042946
6 2 4 0.081655 0.775345 0.090134
7 1 4 0.010043 0.031180 -0.010391
8 0 4 100.000000 NaN 100.000000
0 9 3 NaN 0.530400 0.011832
1 8 3 0.059883 0.473414 0.006330
2 7 3 0.058104 0.479511 0.009540
3 6 3 0.059367 0.539952 0.001254
4 5 3 0.048899 0.599799 -0.008538
5 4 3 0.057369 0.650038 -0.058620
6 3 3 0.052796 0.664637 0.128512
7 2 3 0.055088 0.602155 0.105062
8 1 3 0.052314 0.405939 0.012396
9 0 3 100.000000 NaN 100.000000
0 10 2 NaN 0.291148 -0.022733
1 9 2 0.037852 0.263114 -0.061184
2 8 2 0.031744 0.261888 0.072182
3 7 2 0.026053 0.278429 -0.046420
4 6 2 0.025160 0.308544 0.033588
5 5 2 0.027043 0.330855 0.018518
6 4 2 0.023062 0.355467 0.027696
7 3 2 0.025641 0.383779 -0.142862
8 2 2 0.026757 0.400287 0.036987
9 1 2 0.029170 0.432749 0.042697
10 0 2 0.042096 NaN 6.242737
0 11 1 NaN 0.083634 0.282780
1 10 1 0.011728 0.078903 0.014923
2 9 1 0.010272 0.082378 -0.030417
3 8 1 0.008778 0.089670 -0.017137
4 7 1 0.005230 0.088074 0.001795
5 6 1 0.006741 0.092497 -0.009243
6 5 1 0.009019 0.121387 0.024617
7 4 1 0.007946 0.146089 0.068148
8 3 1 0.008877 0.153441 0.008456
9 2 1 0.007930 0.230180 0.009168
10 1 1 0.010857 0.310193 -0.000974
11 0 1 100.000000 NaN 100.000000
0 12 0 NaN 0.000000 NaN
1 11 0 0.000024 0.000004 NaN
2 10 0 0.000256 0.000484 NaN
3 9 0 0.000146 0.000570 NaN
4 8 0 0.000140 0.000718 NaN
5 7 0 0.000100 0.000520 NaN
6 6 0 0.000116 0.001899 NaN
7 5 0 0.000174 0.004456 NaN
8 4 0 0.000285 0.006070 NaN
9 3 0 0.000171 0.006831 NaN
10 2 0 0.000190 0.012930 NaN
11 1 0 0.000242 0.023072 NaN
12 0 0 0.006223 NaN NaN

Note that in some cases when we have no fundamentalists, the price explodes exponentially, and profits can become infinite or very large. We cap the profits at 100 to deal with these outliers.

We can vary the utility function. Here we consider a simple risk-averse utility function of the form:

RiskAverseUtility
## function(payoff.mean, payoff.variance, coefficient = 2) {
##   payoff.mean - coefficient * payoff.variance
## }

This gives us very different payoffs, which we can see by examining the first six rows of the payoff matrix:

game.abm.ra <- HeuristicGame('data/payoff-matrix.csv', strategies, 1, Utility=RiskAverseUtility)
kable(head(game.abm.ra$payoffs.normalised), col.names=columns)
\(N_C\) \(N_F\) \(N_N\) \(u(e_F)\) \(u(e_C)\) \(u(e_N)\)
0 0 12 NaN NaN -3.270642
0 1 11 NaN -3.509895 -0.157542
1 0 11 -0.382272 NaN 24.390954
0 2 10 NaN 4.361476 -0.043086
1 1 10 0.773708 -0.711277 -0.041043
2 0 10 3.781281 NaN 12.184021

Large-population dynamics

We can now analyse evolution between these strategies, assuming that there is a very large population of brokers who repeatedly trade in randomly-chosen groups of 12 at a time.

For the risk-neutral (rn) case, we obtain:

if (!var.exists('game.abm.rn.analysed'))
  game.abm.rn.analysed <-  Analyse(game.abm.rn, initial.values = initial.values.random, time=seq(0,100,by=0.05), parallel=T)
result <- plot(game.abm.rn.analysed)

For the risk-averse (ra) case, we obtain:

if (!var.exists('game.abm.ra.analysed')) 
  game.abm.ra.analysed <- Analyse(game.abm.ra, initial.values = initial.values.random, time=seq(0,100,by=0.05), parallel=T)
result <- plot(game.abm.ra.analysed)

We can also plot a single trajectory as a time-series:

trajectory <- game.abm.ra.analysed$trajectories[[1]]
plot.ts(trajectory)

Small-population dynamics

Next we will consider co-evolution in a small population with the original 12 agents, and with a Boltzmann selection function: \(p_i = \frac{e^{f_i}/kT}{\sum_i e^{f_i/kT}}\). Initially we consider a temperature of \(T = 1\)

game.abm.ra.ga <- GeneticAlgorithmGame(game.abm.ra, temperature = 1)
if (!var.exists('game.abm.ra.ga.analysed'))
 game.abm.ra.ga.analysed <- Analyse(game.abm.ra.ga, initial.values = initial.values.random,
                                    times = seq(0, 200), parallel=T)
result <- plot(game.abm.ra.ga.analysed)

As we can see, this looks very random. Let’s look at an individual trajectory as a time-series:

trajectory <- game.abm.ra.ga.analysed$trajectories[[1]]
plot.ts(trajectory)

Notice that we see a punctuated equilibrium with high-levels of fundamentalism.

Let’s see what happens when we increase or decrease the temperature.

game.abm.ra.ga.high.temp <- game.abm.ra.ga
game.abm.ra.ga.high.temp$temperature <- 5
if (!var.exists('game.abm.ra.ga.high.temp.analysed'))
  game.abm.ra.ga.high.temp.analysed   <- Analyse(game.abm.ra.ga.high.temp, initial.values = initial.values.random,
                                                  times = seq(0, 200), parallel=T)
result <- plot(game.abm.ra.ga.high.temp.analysed)

Paradoxically, increasing the selection intensity results in the elimination of fundamentalists. Initially the fundamentalists increase in frequency, but then immediately get wiped out:

trajectory <- game.abm.ra.ga.high.temp.analysed$trajectories[[1]]
plot.ts(trajectory)

We can see why this happens by examining the payoffs when either chartists, fundamentalists or noise-traders alternatively dominate the population:

Payoffs to a single agent in a population of chartists:

payoff.against.C <- IndividualPayoff(game.abm.ra.ga.high.temp, c(1, 0, 0))
payoff.against.C
## [1]   0.000937   0.021354 100.000000

Payoffs to a single agent in a population of fundamentalists:

payoff.against.F <- IndividualPayoff(game.abm.ra.ga.high.temp, c(0, 1, 0))
payoff.against.F
## [1] -0.000005  0.000000  0.233001

Payoffs to a single agent in a population of fundamentalists:

payoff.against.N <- IndividualPayoff(game.abm.ra.ga.high.temp, c(0, 0, 1))
payoff.against.N
## [1] -0.382272 -3.509895 -3.270642

This contrasts with the payoffs when there is a more even mix of strategies in the population:

payoff.against.uniform <- IndividualPayoff(game.abm.ra.ga, c(1/3, 1/3, 1/3))
payoff.against.uniform
## [1]  0.7602154  0.8298737 -0.7621616

Comparing in a table:

kable(t(cbind(payoff.against.C, payoff.against.F, payoff.against.N, payoff.against.uniform)), col.names = strategies)
C F N
payoff.against.C 0.0009370 0.0213540 100.0000000
payoff.against.F -0.0000050 0.0000000 0.2330010
payoff.against.N -0.3822720 -3.5098950 -3.2706420
payoff.against.uniform 0.7602154 0.8298737 -0.7621616