1. Our first agent-based evolutionary model

1.4. Analysis of these models

1. Two complementary approaches

Agent-based models are usually analyzed using computer simulation and/or mathematical analysis.

  • The computer simulation approach consists in a) running many simulations –i.e. sampling the model many times– and, with the data thus obtained, b) trying to infer general patterns and properties of the model. In terms of types of reasoning, this approach involves logical deduction first (when the data is generated by the computer following the algorithmic rules that define the model) and induction later (when general patterns are inferred from the generated data).
  • Mathematical approaches do not look at individual simulation runs, but instead analyze the rules that define the model directly, and try to derive their logical implications. Mathematical approaches use deductive reasoning only, so their conclusions follow with logical necessity from the assumptions embedded in the model (and in the mathematics employed).

These two approaches are by no means incompatible; they are complementary in the sense that they can provide fundamentally different insights on the same model and, furthermore, there are plenty of synergies to be exploited by using the two approaches together (Izquierdo et al., 2013).

Here we provide several references to material that is helpful to analyze the agent-based models we have developed in this chapter of the book, and illustrate its usefulness with a few examples. Section 2 below deals with the computer simulation approach, while section 3 addresses the mathematical analysis approach.

2. Computer simulation approach

The task of running many simulation runs –with the same or different combinations of parameter values– is greatly facilitated by a tool named BehaviorSpace, which is included within NetLogo and is very well documented at NetLogo website. Here we provide an illustration of how to use it.

Consider a coordination game defined by payoffs [[3 0][0 2]] and played by 1000 agents who revise their strategy with probability 0.01 in every tick (following the imitate-if-better rule without noise). This model is stochastic and we wish to study how it usually behaves, departing from a situation where both strategies are equally represented. To that end, we could run several simulation runs (say 1000) and plot the average fraction of 0-strategists in every tick, together with the minimum and the maximum values observed across runs in every tick. An illustration of this type of graph is shown in figure 1.

Figure 1. Average proportion of 0-strategists in an experiment of 1000 simulation runs. Orange error bars show the minimum and maximum values observed across the 1000 runs. Payoffs: [[3 0][0 2]]; prob-revision: 0.01; noise 0; initial conditions [500 500].

To set up the computational experiment that will produce the data required to draw figure 1, we just have to go to Tools (in the upper menu of NetLogo) and then click on BehaviorSpace. The new experiment can be set up as shown in figure 2.

Figure 2. Experiment setup in BehaviorSpace

In this particular experiment, we are not changing the value of any parameter, but doing so is straightforward. For instance, if we wanted to run simulations with different values of prob-revision –say 0.01, 0.05 and 0.1–, we would just write:

[ "prob-revision" 0.01 0.05 0.1 ] 

If, in addition, we would like to explore the values of noise 0, 0.01, 0.02 … 0.1, we could use the syntax for loops, [start increment end], as follows:

[ "noise" [0 0.01 0.1] ] ;; note the additional brackets

If we made the two changes described above, then the new computational experiment would comprise 33000 runs, since NetLogo would run 1000 simulations for each combination of parameter values (i.e. 3 × 11).

The original experiment shown in figure 2, which consists of 1000 simulation runs only, takes a couple of minutes to run. Once it is completed, we obtain a .csv file with all the requested data, i.e. the fraction of 0-strategists in every tick for each of the 1000 simulation runs – a total of 1001000 data points. Then, with the help of a pivot table (within e.g. an Excel spreadsheet), it is easy to plot the graph shown in figure 1. A similar graph that can be easily plotted is one that shows the standard error of the average computed in every tick (see figure 3).[1]

Figure 3. Average proportion of 0-strategists in an experiment of 1000 simulation runs. Orange error bars show the standard error. Payoffs: [[3 0][0 2]]; prob-revision: 0.01; noise 0; initial conditions [500 500].

3. Mathematical analysis approach. Markov chains

From a mathematical point of view, agent-based models can be usefully seen as time-homogeneous Markov chains (see Gintis (2013) and Izquierdo et al. (2009) for several examples). Doing so can make evident many features of the model that are not apparent before formalizing the model as a Markov chain. Thus, our first recommendation is to learn the basics of this theory. Karr (1990), Kulkarni (1995, chapters 2-4), Norris (1997), Kulkarni (1999, chapter 5), and Janssen and Manca (2006, chapter 3) are all excellent introductions to the topic.

The models developed in this chapter can be seen as time-homogeneous Markov chains on the finite space of possible strategy distributions. This means that the number of agents that are following each possible strategy is all the information we need to know about the present –and the past– of the stochastic process in order to be able to –probabilistically– predict its future as accurately as it is possible. Thus, the number of possible states in these models is \binom{N + s - 1}{s-1}, where N is the number of agents and s the number of strategies.[2]

In some simple cases, a full Markov analysis can be conducted by deriving the transition matrix of the Markov chain and operating with it. The example in section 2 above is a case in point. In that coordination game, the number of possible states is just  \binom{N + s - 1}{s-1} = \binom{1000 + 2 - 1}{2-1} = 1001, and transition probabilities are not hard to compute (with the help of a computer).

However, most often a full Markov analysis is unfeasible due to one or more of the following –somewhat related– reasons:

  • The number of states is too large.[3]
  • The transition probabilities are hard to derive.
  • The model has too many parameters.

In these common cases where we cannot operate directly with the transition matrix of the Markov chain, can we still say something about the model using mathematical approaches? Certainly so! Even when a full analysis of the Markov chain is not viable, there are often many insights to be gained by focusing on specific aspects of the model. Below we summarize some of the mathematical approaches that can often be conducted.

3.1. Partition of the state space into its closed communicating classes (and all the rest)

In many cases where the transition matrix cannot be easily derived or it is unfeasible to operate with it, one can still partition the state space as the union of all closed communicating classes plus another class containing all the states that belong to non-closed communicating classes.[4] This partition can uncover many properties of the Markov chain even without knowing the exact values of its transition matrix, and these properties can yield neat insights about the dynamics of the associated model. You can find examples of application of this technique to several well-known agent-based models in the social simulation literature in Izquierdo et al. (2009).

As an illustration, consider the coordination game example in section 2. The partition of the state space of this model without noise reveals that simulations will necessarily end up in one of two possible absorbing states: the state where everyone chooses strategy 0, or the state where everyone chooses strategy 1. The probability of reaching one absorbing state or the other depends on initial conditions and can be numerically approximated to any degree of accuracy. In stark contrast, if there is noise, then it is possible to go from any state to any other state, so the whole state space constitutes one single communicating class (and we say that the Markov chain is irreducible). This implies that the probability of finding the system in each of its states in the long run is strictly positive and independent of the initial conditions. It also implies that the limiting distribution coincides with the occupancy distribution. Again, this distribution can be numerically approximated to any degree of accuracy.

3.2. Deterministic approximations of transient dynamics when the population is large. The mean dynamic

When the number of agents is sufficiently large, the mean dynamic of the process provides a good deterministic approximation to the dynamics of the process over finite time spans. The mean dynamic of a Markov chain is derived by taking the limit when the population size goes to infinity, so it characterizes the dynamics of the model over moderate time spans, as long as the population is large enough. For an intuitive introduction to the mean dynamics, see Izquierdo and Izquierdo (2013), who also provide a detailed analysis of the Hawk-Dove game with the imitate-if-better rule. For a concise and rigorous presentation of the most relevant results on Markov chains and on their convergence to the mean dynamics in the context of evolutionary game theory, see Sandholm (2010, chapter 10).

As an illustration of the usefulness of the mean dynamic to approximate transient dynamics, consider the coordination game example in section 2. The mean dynamic for this model reads:

\dot{x} = x (x-1)\left(x^2-3 x+1\right)

where x stands for the fraction of 0-strategists.[5][6] The solution of this ordinary differential equation with initial condition x_0 = 0.5 is shown in figure 4 below. It is clear that the mean dynamic provides a remarkably good approximation to the average transient dynamics plotted in figures 1 and 3.

Figure 4. Trajectory of the mean dynamic of the example in section 2, showing the proportion of 0-strategists as a function of time (rescaled to match figures 1 and 3).

Naturally, the mean dynamic can be solved for many different initial conditions, providing an overall picture of the transient dynamics of the model when the population is large. Figure 5 below shows an illustration, created with the following Mathematica® code:

Plot[
 Evaluate[
  Table[
   NDSolveValue[{x'[t] == x[t] (x[t] - 1) (x[t]^2 - 3 x[t] + 1), 
      x[0] == x0}, x, {t, 0, 10}][t/100]
  , {x0, 0, 1, 0.01}]
 ], {t, 0, 1000}]
Figure 5. Trajectories of the mean dynamic of the example in section 2, showing the proportion of 0-strategists as a function of time (rescaled to match figures 1 and 3) for different initial conditions.

To compare agent-based simulations of the imitate-if-better rule and its mean dynamics in 2×2 symmetric games, you may want to play with the purpose-built demonstration titled Expected Dynamics of an Imitation Model in 2×2 Symmetric Games. To solve the mean dynamic of the imitate-if-better rule in 3-strategy games, you may want to use this demonstration.

3.3. Diffusion approximations to characterize dynamics around equilibria

“Equilibria” in finite population dynamics are often defined as states where the expected motion of the (stochastic) process is zero. Formally, these equilibria correspond to the rest points of the mean dynamic of the original stochastic process. Thus, at such equilibria, the expected flows of agents switching between different strategies cancel one another out but, naturally, agents keep revising and changing strategies, potentially in a stochastic fashion.

As an example, consider a Hawk-Dove game with payoffs [[0 3][1 2]] and the imitate-if-better revision rule without noise. The mean dynamic of this model is:

\dot{x} = x (1-x)(1-2 x)

where x stands for the fraction of 0-strategists, i.e. “Hawk” agents.[7] Solving the mean dynamic reveals that most large-population simulations starting with at least one “Hawk” and at least one “Dove” will tend to approach the state where half the population play “Hawk” and the other play “Dove”, and stay around there for long. Figure 6 below shows several trajectories for different initial conditions.

Figure 6. Trajectories of the mean dynamic of an imitate-if-better Hawk-Dove game, showing the proportion of “Hawks” as a function of time for different initial conditions.

Naturally, simulations do not get stuck in the half-and-half state, since agents keep revising their strategy in a stochastic fashion  (see figure 7). To understand this stochastic flow of agents between strategies near equilibria, it is necessary to go beyond the mean dynamic. Sandholm (2003) shows that –under rather general conditions– stochastic finite-population dynamics near rest points can be approximated by a diffusion process, as long as the population size N is large enough. He also shows that the standard deviations of the limit distribution are of order \frac{1}{\sqrt{N}}.

To illustrate this order \frac{1}{\sqrt{N}}, we set up one simulation run starting with 10 agents playing “Hawk” and 10 agents playing “Dove”. This state constitutes a so-called “Equilibrium”, since the expected change in the strategy distribution is zero. However, the stochasticity in the revision protocol and in the matching process imply that the strategy distribution is in perpetual change. In the simulation shown in figure 7, we modify the number of players at runtime. At tick 10000, we increase the number of players by a factor of 10 up to 200 and, after 10000 more ticks, we set n-of-players to 2000 (i.e., a factor of 10, again). The standard deviation of the fraction of players using strategy “Hawk” (or “Dove”) during each of the three stages in our simulation run was: 0.1082, 0.0444 and 0.01167 respectively. As expected, these numbers are related by a factor of approximately \frac{1}{\sqrt{10}}.[8]

Figure 7. A simulation run of an imitate-if-better Hawk-Dove game, set up with 20 agents during the first 10000 ticks, then 200 agents during the following 10000 ticks, and finally 2000 agents during the last 10000 ticks. Payoffs: [[0 3][1 2]]; prob-revision: 0.01; noise 0; initial conditions [10 10].

3.4. Stochastic stability analyses

In the last model we have implemented in this unit, recall that when noise is strictly positive, it is possible to go from any state to any other state, so the model’s infinite-horizon behavior is characterized by a unique stationary distribution regardless of initial conditions (see section 3.1 above). This distribution has full support (i.e. all states will be visited infinitely often) but, naturally, the system will spend much longer in certain areas of the state space than in others. If the noise is sufficiently small (but strictly positive), the infinite-horizon dynamics of the Markov chain are often concentrated in just a few states. Stochastic stability analyses are devoted to identifying such states, which are often called stochastically stable states, and are a subset of the absorbing states of the process without noise.

To learn about this type of analysis, the following references are particularly useful: Vega-Redondo (2003, section 12.6)Fudenberg and Imhof (2008), Sandholm (2010, chapters 11 and 12) and Wallace and Young (2015).

To illustrate the applicability of stochastic stability analyses, consider our imitate-if-better model where agents play a Prisoner’s Dilemma with payoffs [[2 4][1 3]] and strictly positive noise. It can be proved that the only stochastically stable state in this model is the state where everyone defects.[9] This means that, as noise tends to 0, the infinite-horizon dynamics of the model will concentrate on that single state.

4. Exercises

Exercise 1. Derive the mean dynamic of a Prisoner’s Dilemma game for the imitate-if-better protocol.

Exercise 2. Derive the mean dynamic of the coordination game discussed in section 0.1 (with payoffs [[1 0][0 2]]) for the imitative pairwise-difference protocol.

Exercise 3. Derive the mean dynamic of the coordination game discussed in section 0.1 (with payoffs [[1 0][0 2]]) for the best experienced payoff protocol.


  1. The standard error of the average equals the standard deviation of the sample divided by the square root of the sample size. In our example, the maximum standard error was well below 0.01.
  2. This result can be easily derived using a "stars and bars" analogy.
  3. As an example, in a 4-strategy game with 1000 players, the number of possible states in our model is \binom{1000 + 4 - 1}{4-1} = 167668501.
  4. This partition is unique (see decomposition theorem by Chung (1960)).
  5. For details, see Izquierdo and Izquierdo (2013).
  6. One unit of clock time is defined in such a way that each player expects to receive one revision opportunity per unit of clock time. In this particular example, since prob-revision = 0.01, one unit of clock time corresponds to 100 ticks (i.e. 1 / prob-revision).
  7. For details, see Izquierdo and Izquierdo (2013).
  8. A good approximation for these standard deviations over long time spans when the number of agents is large can be exactly computed (see Izquierdo et al. (2018, example 3.1)).
  9. The proof can be conducted using the concepts and theorems put forward by Ellison (2000). Note that the radius of the state where everyone defects is 2 (i.e. 2 mutations are needed to leave its basin of attraction), while its coradius is just 1 (one mutation is enough to go from the state where everyone cooperates to the state where everyone defects).

License

Icon for the Creative Commons Attribution 4.0 International License

Agent-Based Evolutionary Game Dynamics by Luis R. Izquierdo, Segismundo S. Izquierdo & William H. Sandholm is licensed under a Creative Commons Attribution 4.0 International License, except where otherwise noted.

Share This Book