## Abstract

In this paper, we present different applications of finite state mean field games to socio-economic sciences. Examples include paradigm shifts in the scientific community or consumer choice behaviour in the free market. The corresponding finite state mean field game models are hyperbolic systems of partial differential equations, for which we present and validate different numerical methods. We illustrate the behaviour of solutions with various numerical experiments, which show interesting phenomena such as shock formation. Hence, we conclude with an investigation of the shock structure in the case of two-state problems.

## 1. Introduction

Mean field games have become a powerful mathematical tool to model the dynamics of agents in economics, finance and the social sciences. Different settings have been considered in the literature such as discrete and continuous in time or finite and continuous state space. Originally, finite state mean field games [1–6] were studied as an attempt to understand the more general continuous state problems introduced by Lasry & Lions [7–9] as well as by Huang *et al.* [10,11]. For additional information, see also the recent surveys such as [12–16].

In this paper, we apply finite state mean field games to two classical problems in socio-economic sciences: consumer choice behaviour and paradigm shifts in a scientific community [17]. These mean field game models give rise to systems of hyperbolic partial differential equations, which serve as a starting point for the proposed numerical methods.

The mathematical modelling is based on the following situation. Let us consider a system of *N*+1 identical players or agents, which can switch between different states. Each player is in a state and can choose a switching strategy to any other state . The only information available to each player, in addition to its own state, is the number *n*_{j} of players it sees in the different states . The fraction of players in each state is denoted by *θ*_{i}=*n*_{i}/*N*, and we define the probability vector *θ*=(*θ*_{1},…,*θ*_{d}), , which encodes the statistical information about the ensemble of players. Each player faces an optimization problem over all possible switching strategies. Here, the key question is the existence of a Nash equilibrium and to determine the limit as the number of players tends to infinity. Gomes *et al.* [5] showed that the *N*+1 player Nash equilibrium always exists. The Nash equilibrium in the limit (called mean field Nash equilibrium) satisfies certain ordinary differential equations, which may not have unique solutions. Because of this non-uniqueness, the development of numerical methods, that capture the appropriate equilibrium, is of high scientific interest. We recall that the value function *U*^{i}=*U*^{i}(*θ*,*t*), for a player in state *i* in the limit , satisfies the hyperbolic system
1.1
at least for short time (see Gomes *et al.* [5]). Note that the distribution of players in the different states is given by *θ*, , , , and , where ∂_{θj} denotes the partial derivative with respect to the variable *θ*_{j}.

In general, first-order hyperbolic equations do not admit smooth solutions, and one needs to consider an appropriate notion of such. Adequate definitions of solutions are well known for conservation laws and equations which admit a maximum principle, e.g. Hamilton–Jacobi equations. Until now, the appropriate notion of solutions for (1.1), which encodes the mean field limit, has not been clear.

In this paper, we present a numerical method, which is based on the Nash equilibrium equations for *N*+1 agents. Therefore, these equations will, in the case of convergence, automatically yield the appropriate limit. Convergence for short times was shown in [5]. For certain classes of finite state mean field games, called potential mean field games, system (1.1) can be regarded as the gradient of a Hamilton–Jacobi equation [1,5]. We use this property to validate the presented numerical method. Our computational experiments show the expected formation of shocks. We analyse these shock structures by introducing an auxiliary conservation law. This allows us to derive a Rankine–Hugoniot condition that characterizes the qualitative behaviour of such systems.

This paper is organized as follows: we start §2 by recalling the *N*+1 player model, which gives rise both to (1.1) and to the numerical method presented here. Section 3 focuses on two-state mean field games, the numerical implementation of (1.1) and our applications in socio-economic sciences. We illustrate the behaviour of the various models in §4. In §5, we briefly analyse the shock structure for two-state mean field games.

## 2. Finite state mean field games

We start with a more detailed presentation of finite state mean field games and the formal derivation of (1.1). The model presented here extends slightly the one developed in [5] as it allows for interactions between agents as described in the next section. Furthermore, the *N*+1 player problem gives rise to a discretization of (1.1). This discretization is the basis of our numerical method detailed in §4.

### (a) *N*+1 player problem

We consider a system of *N*+1 identical players or agents. Let us fix one of them, referred to as the *reference player*, and denote by **i**^{t} its state at time *t*. All other players can be in any state at each time *t*. We denote by the number of players (distinct from the reference player) at state *j* at time *t*, and . Players change their states according to two mechanisms: either by a Markovian switching rate chosen by the player, or by interactions with the other players. The interactions are determined by the following simplified model. We assume that any pair of players will meet at random times whose distribution is exponential. As an outcome of the meeting, one of the two players chooses to join the state of the other player. This is a typical behaviour, for instance, in opinion dynamics, where players change their opinions by interactions with other players. We assume for the moment that all players (except the reference player) have chosen an identical Markovian strategy . The reference player is allowed to choose a possibly different strategy . Given these strategies, the joint process (**i**^{t},**n**^{t}) is a time inhomogeneous Markov chain (neither **n**^{t} nor **i**^{t} is Markovian when considered separately). It is sufficient to specify this Markov process, as in the classical case of Markov chains with discrete space state and discrete time, by defining the generator of the process. We split the generator into three parts, each one corresponding to different interactions, which will be discussed later on. Hence, we have

Let *e*_{k} be the *k*th vector of the canonical basis of and *e*_{jk}=*e*_{j}−*e*_{k}. Then
The terms and correspond to the transitions owing to the switching strategies *α* and *β*. The first two terms of the generator are similar to those described in more detail in [5], but we briefly recall their interpretation for completeness.

Select one of the players distinct from the reference player. Denote by **k**^{t} its position at time *t*, and by the process that records the number of other players in any state from the point of view of this player. Suppose further that there are no interactions (). Then, for *j*≠*k*, we have for any time increment *δ* small enough, that
Assuming symmetry and independence of transitions from any state *k* to a state *j*, *k*≠*j*, we have
where the transition rates of the process **n**^{t} are given by
2.1
Similarly, the reference player switching probabilities are

The transitions between different states owing to interactions give rise to the term . Its particular structure comes from the assumption that any two players with distinct states *j* and *k* can meet with a rate *ω*_{jk}/*N* (with *ω*_{kj}=*ω*_{jk}≥0). As a result of this interaction, both end in either state *j* or state *k* (with probability respectively).

We assume that all players have the same running cost determined by a function as well as an identical terminal cost *U*_{T}(*θ*), which is Lipschitz continuous in *θ*. The running cost *c*(*i*,*θ*,*α*) depends on the state *i* of the player, the mean field *θ*, that is the distribution of players among states, and on the switching rate *α*. As in [5], we suppose that *c* is Lipschitz continuous in *θ* with a Lipschitz constant (with respect to *θ*) bounded independently of *α*. Let the running cost *c* be differentiable with respect to *α*, and (∂*c*/∂*α*)(*i*,*θ*,*α*) be Lipschitz with respect to *θ*, uniformly in *α*. We assume that for each , the running cost *c*(*i*,*θ*,*α*) does not depend on the *i*th coordinate *α*_{i} of *α*. Furthermore, we make the following additional assumptions on the cost *c*:

(A1) For any , , , with

*α*_{j}≠*α*_{j}′, for some*j*≠*i*,(A2) The function

*c*is superlinear in*α*_{j},*j*≠*i*, i.e.

Let us fix a reference player and set the Markovian strategy *β* for the remaining *N* players. The objective of the reference player is to minimize its total cost. The minimum over all Markovian strategies *α* is given by
2.2

Define for , Δ_{i}*φ*:=(*φ*^{1}−*φ*^{i},…,*φ*^{d}−*φ*^{i}) and the generalized Legendre transform of *c* by
2.3
Note that *h* depends only on the differences between the coordinates of the variable *z*, that is, if then .

The function is the solution to the ODE
Next, we define, for *j*≠*i*,
2.4
If *h* is differentiable, for *j*≠*i*,
2.5
For convenience and consistency with (2.5), we require that
2.6
Then, the optimal strategy for the reference player is given by

We say that a strategy *β* is a Nash equilibrium if the optimal response of the reference player is *β* itself, i.e. . Thus, by setting , we obtain the Nash equilibrium equation for the value function, i.e.
2.7
where we define
2.8

### (b) Formal asymptotic behaviour

Next, we investigate the asymptotic behaviour of the *N*+1 player dynamics (2.7) as the number of players tends to infinity. Hence, we suppose that there is a smooth function such that

Then, we have the following expansions:

We observe that for fixed *j* and *k* the operators and *ω*_{jk} *θ*_{j}*θ*_{k}(∂_{θj}−∂_{θk})^{2} are degenerate elliptic operators. The first one is degenerate, because *θ*_{k}, ; the second, because *ω*_{jk}≥0. Hence, their sum is also a degenerate operator. Therefore, the combination of the second-order terms in the expansion of and can be written as
for a suitable non-negative matrix *b*. We conclude that (2.7) can be formally approximated by the parabolic system
2.9
with suitable . Note that the viscosity term on the right-hand side of (2.9) does not come from an artificial argument in order to solve (2.9) numerically or to obtain better regularity results. On the contrary, this term is built in the original *N* player model. Furthermore, *g*^{N} converges locally uniformly in compacts to
2.10

Therefore, the limit of (2.9) is (1.1). Note that (1.1) does not depend on the interactions *ω*_{jk} between players. Furthermore,
2.11
because , from (2.6). Additionally,
2.12
using (2.4) in equation (2.10).

### (c) Potential mean field games

Next, we consider a special class of mean field games, in which system (1.1) can be written as the gradient of a Hamilton–Jacobi equation. Suppose that
2.13
and *f*(*i*,*θ*)=∂_{θi}*F*(*θ*), for some potential . We set
2.14

Let be a continuous function and consider a sufficiently smooth solution to the Hamilton–Jacobi equation
2.15
Note that . In some cases, it is possible to reduce the dimensionality at the price of introducing suitable boundary conditions (see §4*b*). This reduction is used in the applications presented later.

Set
2.16
If we differentiate (2.15) with respect to *θ*_{i}, we obtain
The first term on the right-hand side can be written as
taking into account the identity ∂_{θi}*U*^{j}=∂_{θj}*U*^{i}. From this, we obtain that
and deduce, using (2.13), that *U*^{i} is indeed a solution of (1.1).

Finally, we remark that potential mean field games have similar properties and close connections to calculus of variations. For instance, long-time convergence properties of these problems can be addressed through *Γ*-convergence techniques (e.g. [6]).

### (d) Characteristics and shocks

We conclude this section by discussing the characteristic ODE for (1.1)
2.17
for , *t*∈[0,*T*], and . As presented in [5], the solutions to this ODE satisfy the initial–terminal boundary conditions
2.18
and correspond precisely to the mean field Nash equilibria for the limit for an initial distribution of players given by . In general, system (2.17)–(2.18) admits multiple solutions, hence the Nash equilibrium is not unique. Multiple solutions to the characteristic ODE correspond to the existence of shocks in (1.1). The importance of understanding shock formation and the corresponding relevant solutions of (1.1), which encode the mean field limit, cannot be underestimated and it is essential to select the appropriate solutions to (2.17)–(2.18).

## 3. Two-state mean field games in socio-economic sciences

Here, we present two model applications of finite state mean field games in socio-economic sciences. For reasons of clarity and readability, we focus only on two-state problems. With minor modifications, our approach can be adapted to a larger number of states. We start by stating the explicit equations for the case where agents can choose between two options. Next, we address two explicit examples, which illustrate the strength of finite state mean field games methods in socio-economic modelling: first, consumer choice behaviour; second, a paradigm shift model in the scientific community. These examples shall serve as a starting point for future research in this direction.

### (a) Two-state problems

Consider a two-state mean field game, where the fraction of players in state 1 or state 2 is given by *θ*_{i}, *i*=1,2 with *θ*_{1}+*θ*_{2}=1, and *θ*_{i}≥0. Because the limit equation (1.1) does not depend on the interactions (although the *N*+1 player model does), we set *ω*=0. Note that *ω*≠0 would result in different numerical methods (and potentially different solutions) for (1.1). We suppose further that the running cost *c*=*c*(*i*,*θ*,*μ*) in (2.2) depends quadratically on the switching rate *μ*, i.e.
3.1

Then, 3.2

The optimal switching rate *α** is given by

Because we conclude from (2.10) that 3.3

Note that if the function *f* is a gradient field, i.e. *f*=∇*F*, the two-state problem is a potential mean field game, cf. §2*c*. In this case, for , the function *H*, defined by (2.14), reads as
3.4

The above calculations allow us to introduce a numerical method based on the two-state mean field model for *N*+1 players. Let *n*_{i}, *i*=1,2 denote the number of players in state *i* (as seen by the reference player excluding itself) and *N*=*n*_{1}+*n*_{2}. The vector *n* gives the number of players in each state, i.e. *n*=(*n*_{1},*n*_{2})=(*n*_{1},*N*−*n*_{1}). As in (2.8), we have
Then, equation (2.7) for the value function reads as
3.5
which can be rewritten as
3.6

Because *θ*_{1}+*θ*_{2}=1, we introduce the new variable *ζ*∈[0,1]; hence, *θ*=(*ζ*,1−*ζ*). We split the domain [0,1] into *N* equidistant subintervals and define *ζ*_{k}=*k*/*N*, . The variable *ζ*_{k} corresponds to the fraction of players in state 1. Then, the fraction of players in state 2 is given by 1−*ζ*_{k}=(*N*−*k*)/*N*. Consequently, (3.6) takes the form
3.7

Note that in system (3.7) the terms *ζ*_{k} and (1−*ζ*_{k}) vanish for *k*=0 and *k*=*N*, thus no particular care has to be taken concerning the ghost points at *ζ*_{N+1} and *ζ*_{−1}. This is the discrete analogue to not imposing boundary conditions on (1.1). A similar situation occurs in state-constrained problems for Hamilton–Jacobi equations.

### (b) Paradigm shift

According to Kuhn [18], a paradigm shift corresponds to a change in a basic assumption within the ruling theory of science. Classical cases of paradigm shifts are the transition from Ptolemaic cosmology to Copernican one, the development of quantum mechanics which replaced classical mechanics on the microscopic scale or the acceptance of Mendelian inheritance as opposed to pangenesis.

A common assumption in theoretical models for research dynamics is the fact that scientists are rewarded from recognition from others [19,20]. This recognition may take the form of citations, prizes or other financial incentives and is further amplified by scientific activity, such as conferences and collaborations, in this field. Bensancenot & Dogguy [17] modelled a paradigm shift in a scientific community by a two-state mean field game approach and analysed the competition between two different scientific hypotheses. In our example, we consider a simpler model, but follow their general ideas and assumptions.

Let us consider a scientific community with *N* researchers working on two different hypotheses. Each researcher working on paradigm *i*, *i*=1,2, wants to maximize his/her productivity measured by a cost function of the form (3.1). Here, the function *f*=*f*(*i*,*θ*) corresponds to the productivity of a researcher working on paradigm *i*, and to the cost of switching to the other objective. Note the negative sign of the switching costs, because agents want to maximize their productivity. We assume that the productivity is directly related to the number of researchers working on the paradigm, because, for example, more scientific activities such as conferences and collaborations take place. In the case of two different fields, *θ*_{1} gives the fraction of researchers working on paradigm 1 and *θ*_{2}=1−*θ*_{1} on paradigm 2. We choose the functions *f*(*i*,*θ*), *i*=1,2, of the form
3.8
These functions are called *productivity functions with constant elasticity of substitution* and are commonly used in economics to combine two or more productive inputs (in our case, scientific activities in the different fields) to an output quantity. The constant , denotes the *elasticity of substitution*, and it measures how easy one can substitute one input for the other. The constants *a*_{i}∈[0,1] measure the dependence of paradigm *i* with respect to the other. If *a*_{i} is close to one, the field is more autonomous and little influenced by the activity in the other field. The choice of *f* is also motivated by the assumption that scientists place greater weight on theories which are accepted by a larger community than otherwise.

### (c) Consumer choice

Consumer choice models relate preferences to consumption expenditure. Research has shown that consumer behaviour is a very complex process, which involves elements from psychology, sociology, marketing and economics. We propose a very stylized model, in which consumer preferences are solely based on the price of the goods available and the fraction of consumers using them. This model may be applicable to study the dynamics in the mobile phone sector, where consumer choices are strongly influenced by the number of people using the same provider and the costs of the contract. We also refer to [21] for other economic applications of two-state mean field games.

Let us consider two choices of consumption goods and denote by *θ*_{1} the fraction of agents consuming good 1 and by *θ*_{2}=1−*θ*_{1} the fraction consuming good 2. We assume that the price of a good is strongly determined by the consumption rate; in particular we choose, for *i*=1,2,
3.9
where corresponds to the minimum price of the good. In economic literature, the function *f* is called the *isoelastic utility function*. It expresses the utility in terms of the consumption *θ* and exhibits a constant risk aversion *η*. Risk aversion measures the reluctance of agents to accept a bargain with an uncertain pay-off in place of another bargain with a more certain, but possibly lower, expected pay-off.

## 4. Numerical simulations

In the following, we illustrate the behaviour of the proposed numerical method for different applications. Owing to the hyperbolic nature of (1.1) shocks may arise even in the case of smooth terminal conditions. As discussed in §2*d*, shocks in the utility function are caused by the non-uniqueness of the Nash mean field equilibrium. When shocks are present, a small variation in the initial distribution of players may give rise to a substantially different behaviour of the mean field through the ODEs (2.17)–(2.18). Hence, shock formation in the utility function indicates the existence of a clear threshold for the preference of a state in terms of the distribution of players.

Let *N*=100, i.e. the interval [0,1] is discretized into 100 equidistant intervals. Each grid point corresponds to the percentage of players being in state 1. System (3.7) is solved using an explicit in time discretization with time steps of size Δ*t*=10^{−4}. In all examples in this section the terminal time *T* is set to *T*=10 if not stated otherwise.

### (a) Numerical examples

*Example I (shock formation)*: in this first example, we illustrate the formation of shocks, a phenomenon well known for Hamilton–Jacobi equations. We choose a terminal cost of the form
a running cost as in (3.1) with *f*(1,*θ*)=1−*θ*_{1} and *f*(2,*θ*)=1−*θ*_{2}=*θ*_{1}. Figure 1 clearly illustrates the formation of a shock for smooth terminal data. This shock is also evident when we consider the difference *U*^{1}−*U*^{2} of the utilities. This difference is a relevant variable in this problem, because both *g* and *h*, given by (2.12) and (2.3) respectively, depend only on the difference between the utilities. This structure is explored in more detail in §5.

*Example II (paradigm shift)*: in this example, we illustrate the outcome of a two-state mean field game modelling a paradigm shift (§3*b*) within a scientific community. Note that we use the negative cost functional in this example because we always consider minimization problems. The terminal utilities are given by
4.1
and the parameters in (3.8) are set to , . In figure 2, we observe the paradigm shift within the scientific community. At *T*=10, the optimal states are *θ*_{1}=1 and *θ*_{2}=1, because the functions *U*^{1} and *U*^{2} take their minimum value at these points, respectively. In figure 2, we observe that this is not the case at *t*=0. Here, *U*^{1} takes its minimum value at *θ*_{1}=0, i.e. paradigm 1 is not popular any more.

*Example III (consumer choice)*: in our final example, we consider the consumer choice behaviour (see §3*c*). We set the final utility function to
At time *T*=10, the utility functions take their minimum value at *θ*_{1}=1 and *θ*_{2}=0, i.e. their minimum corresponds to the case that either all of them choose product 1 or product 2, respectively. Figure 3 illustrates the utility functions for two sets of parameters, namely
We observe for the second set of parameters that *U*^{1} takes its minimum value on the interval *θ*_{1}∈[0.65,1]. The jump at *θ*_{1}≈0.65 in both utilities indicates the existence of a critical acceptance rate. If less than 65% of the consumers buy product 1, the price is increasing and the product will not be competitive.

### (b) Potential mean field games

In order to validate our methods, we consider our previous example I that can be written as a potential mean field game. We remark that the numerical method proposed in this paper can be applied to mean field games that do not admit a potential formulation. In example I, we have that *F* and *H* in (3.4) are given by

Then, we compare the numerical simulations of (3.7) with those for the corresponding Hamilton–Jacobi equation as we explain in what follows.

For *i*=1,2, set

In order to simplify the numerical implementation, we perform a dimensionality reduction. Define (*θ*_{1},*θ*_{2})=(*ζ*,1−*ζ*), with *ζ*∈[0,1]. Let *Ψ* be the solution to (2.15) and define
4.2
We observe that
4.3
where is given by

We use Godunov's method and an explicit Runge–Kutta method to discretize (4.3). Particular care has to be taken at the boundary. Because *ζ*∈[0,1] represents the first component of a probability vector, the natural boundary conditions for this problem are state constraints. A possibility to implement this is by supplementing (4.3) with large Dirichlet boundary values, i.e.
4.4

To implement the Dirichlet boundary conditions, we follow the works of Abgrall [22] and Waagan [23]. Again, we consider an equidistant discretization of the interval [0,1] into *N* subintervals of size Δ*ζ*, and we approximate the solution *Υ*(*ζ*,*τ*) to (4.3) by *Υ*^{τ}(*ζ*) for *τ*=*T*−*lΔt*, and *ζ*=*kΔζ*, 0≤*k*≤*N*; . We set *Υ*^{T}(*ζ*)=*Ψ*_{0}(*ζ*,1−*ζ*). Then, the Godunov scheme can be written as
4.5
where *δ*^{+}_{ζ} and *δ*^{−}_{ζ} are the difference operators
and in (4.5) is given by
At the boundary *ζ*=0,1, we set
where and .

If we denote ∂*Υ*/∂*ζ* by ∂*Υ*(*ζ*,*t*)/∂*ζ*=(∂/∂*θ*_{1}−∂/∂*θ*_{2})*Ψ*|_{(ζ,1−ζ)}, we have that
by using (4.2) and (2.16); this provides us a method of evaluating the difference *U*^{1}−*U*^{2} via the potential version of the mean field game. Figure 4 shows the derivative of *Υ* with respect to *ζ*, i.e. the difference *U*^{1}−*U*^{2} evaluated via the potential version, as well as the difference *U*^{1}−*U*^{2} calculated via our numerical method from (3.7) at time *t*=0. The same spatial and temporal discretization (i.e. *N*=100 and Δ*t*=10^{−4}) was used in both simulations.

## 5. Shock structure for two-state problems

Finally, we perform a more detailed investigation of the shock structure in the case of two-state problems. For this purpose, we perform a reduction of the dimension (because the key quantities depend on the difference of the utility functions only) and obtain a hyperbolic scalar equation. Then, we introduce a related conservation law, which yields a Rankine–Hugoniot condition for possible shocks. This new formulation allows us to study finite state mean field games from another numerical perspective and gives new insights into the shock structure and the qualitative behaviour of solutions to (1.1).

### (a) Reduction to a scalar problem

Let *U*(*θ*,*t*) be a *C*^{1} solution to (1.1) with *d*=2. We define *w*(*ζ*,*t*)=*U*^{1}(*θ*,*t*)−*U*^{2}(*θ*,*t*), where *θ*=(*ζ*,1−*ζ*). From (1.1), we have that
5.1

Then, *h* and *g*, given by (2.3) and (2.12), can be written as *h*(*U*^{1},*U*^{2},*θ*_{1},*θ*_{2},*i*)=*h*(*w*(*ζ*,*t*),0,*ζ*,1−*ζ*,*i*) and *g*_{1}(*U*^{1},*U*^{2},*θ*_{1},*θ*_{2})=*g*_{1}(*w*(*ζ*,*t*),0,*ζ*,1−*ζ*). For two-state problems, equation (2.11) gives *g*_{2}=−*g*_{1}. Hence, we obtain

Define *r* and *q* by
and denote ∂*w*/∂*ζ* by ∂*w*(*ζ*,*t*)/∂*ζ*=(∂/∂*θ*_{1}−∂/∂*θ*_{2})(*U*^{1}−*U*^{2})|_{(ζ,1−ζ)}. Then, equation (5.1) for the difference of *U*^{i} can be written as
5.2

Hence, we obtain a single one-dimensional hyperbolic equation for *w*, which describes the evolution of the difference *w*=*U*^{1}(*θ*,*t*)−*U*^{2}(*θ*,*t*). Equation (5.2) serves as a basis for analysing the shock structure in more detail.

### (b) A conservation law and the Rankine–Hugoniot condition

Let us consider the following conservation law associated with (5.2) on the interval [0,1]:
5.3
supplemented with the boundary condition *P*(0,*t*)=*P*(1,*t*)=0 for all times *t*∈[0,*T*].

If *P* is a sufficiently smooth solution to (5.3) and *P*(*ζ*,0)≥0, then the maximum principle implies that *P*(*ζ*,*t*)≥0. Furthermore, if , we have that . Assuming that *P*(*ζ*,0) is a probability distribution, we can regard *P*(*ζ*,*t*) as a probability distribution on the set , , as we have a natural identification for two-state problems: . Hence, equation (5.3) describes an evolution of a probability distribution on . Uncertainty in the initial distribution of the mean field *ζ* can be encoded in the initial condition *P*(*ζ*,0) and propagated through (5.3).

Because equation (5.3) may not have globally smooth solutions, we use the Rankine–Hugoniot condition to characterize certain possibly discontinuous solutions. Let *s*:[0,*T*]→[0,1] be a *C*^{1} curve and suppose that *P* is a *C*^{1} function on both 0<*ζ*<*s*(*t*) and *s*(*t*)<*ζ*<1, for *t*∈[0,*T*]. Assume further that (5.3) holds in this set. Let . We denote by [*B*] the jump of *B* across the curve *s*(*t*), that is, [*B*]=*B*(*s*^{+}(*t*),*t*)−*B*(*s*^{−}(*t*),*t*). Equation (5.3) leads to the Rankine–Hugoniot condition of the form
5.4

If we start with initial condition *P*(*ζ*,0)=1, then the support of *P* is the closure of the set of all mean field states which can be reached from some initial mean field state (all possible choices of *ζ* at time 0). Suppose . Suppose that there is a discontinuity in *P* at the boundary of the set *P*=0. Then, we conclude from the Rankine–Hugoniot condition that

Hence, is a characteristic for (5.2). From the discussion above, it is clear that the solution to (5.3) with initial condition *P*(*ζ*,0)=1 will track the transport of a uniform measure by the characteristics to (5.2). Therefore, numerical simulations of this quantity provide an indicator whether characteristics run into shocks or not.

Using (5.3), we can also derive local Lipschitz bounds for the solution to (5.2):

### Proposition 5.1

*There exist a time t*_{0}<*T and a constant C, depending only the* *norm and the Lipschitz constant of w*(*x*,*T*), *such that for any C*^{1} *solution* *to* (*5.2*), *we have for t*_{0}<*t*≤*T the following estimate*:

### Proof.

Let *P* be a solution to (5.3) with *P*(0,*t*)=*P*(1,*t*)=0, *t*≤*T*, *P*(*ζ*,*t*)≥0, with . We differentiate equation (5.2) with respect to *ζ*, multiply by 2∂_{ζ}*w* and obtain
using that *S*(*ζ*,*t*)=(∂_{ζ}*w*(*ζ*,*t*))^{2}. Then, we can deduce the following estimate:
where we use the fact that *w* is bounded (see [5]) in the last step. Then
Taking *P*(*ζ*,*t*) to be an arbitrary probability measure on [0,1] this yields that
This estimate does not give global bounds, but a nonlinear version of Gronwall's inequality yields the existence of a time *t*_{0}<*T*, such that the estimate holds. □

### (c) Numerical analysis of the shock structure

Finally, we discuss the particular formulation of equations (5.2) and (5.3) as well as the numerical realization of example I presented in §4. Equation (5.2) can be written as
5.5
using that *h* and *g* are given by (3.2) and (3.3), respectively. Similarly, (5.3) takes the form
5.6

Note that the difference *w*=*U*^{1}−*U*^{2} has to satisfy *w*(0,*t*)≥0 and *w*(1,*t*)≤0. The violation of this assertion corresponds to the case of no agents being in state 1, i.e. *ζ*=0 in a situation where state 1 would be preferred to state 2 (and analogously for *ζ*=1). Therefore, equation (5.6) does not require any boundary conditions, because, for *ζ*=0 and *ζ*=1, the advection term vanishes, i.e. ((1−2*ζ*)|*w*|−*w*)/2=0.

We simulate (5.6) using an upwind finite difference scheme with the same parameters as in §4. The initial datum *P*(*ζ*,0) is set to

Owing to the vanishing advection terms at *ζ*=0,1, we set homogeneous Dirichlet boundary conditions for *P*, that is *P*(*ζ*,*t*)=0 if *ζ*∈{0,1}. We choose an equidistant discretization of *N*=125 points and time steps Δ*t*=10^{−4}. The evolution of *P* for example I presented in §4 is illustrated in figure 5. We observe that the function *P* vanishes on a neighbourhood of the shock in *w* (figure 1), which is located at *ζ*=0.5. This suggests that the appropriate notion of weak solution to (5.2) has the property that characteristics can cross once only at the original point. This is similar with what happens for entropy solutions for conservation laws or viscosity solutions of Hamilton–Jacobi equations.

## 6. Conclusion

In this paper, we present effective numerical methods for two-state mean field games and discuss a number of illustrative examples in socio-economic sciences. We compare our method with numerical results obtained from classical and well established schemes for Hamilton–Jacobi equations. If the original equation (1.1) can be written as a potential mean field game, our numerical scheme gives identical results to the ones obtained from discretizing the Hamilton–Jacobi equation itself. This is a quite interesting observation, because there is no reason *a priori* which would imply that the solutions to (1.1) for potential mean field games agree with the gradients of viscosity solutions to (2.15). Indeed, it is conceivable that a different notion of weak solution rather than viscosity solutions could be more appropriate to describe mean field problems. In the case where the original equation does not admit a potential, our method provides a new way to approximate the solution to (1.1). Furthermore, this method is an effective way to select among all possible solutions to the mean field equilibrium equations (2.17)–(2.18) those that arise in the limit . In the general case, when the mean field game does not have a potential structure, it would be important to compare our numerical method with standard schemes for Hamilton–Jacobi equations or conservation laws, such as those in [24,25]. Of course, owing to the special structure of mean field games, these schemes would have to be suitably modified as they do not apply directly. We plan to address these issues in a forthcoming publication.

As presented in the examples, our method captures shocks effectively. We analyse the shock structure using an associated conservation law and prove a local Lipschitz estimate for the solutions of (1.1). Furthermore, we have shown that the proposed numerical method can be used even in the case of non-uniqueness of mean field equilibria, because it is automatically compatible (by its own nature) with the limit . This is, to the best of our knowledge, the only method available in the literature to address problems without uniqueness of mean field Nash equilibrium and to provide a selection criterion for (2.17)–(2.18). Shock formation is not only a challenging issue in mathematical terms, but it also gives interesting insights into the dynamic behaviour of interacting agents. It corresponds to rapid changes in the utility function and indicates the existence of thresholds in the distribution of agents, where the market behaviour changes significantly. In addition, its presence indicates the lack of uniqueness of a mean field Nash equilibrium.

Several challenging and interesting mathematical questions will be addressed in a future paper, such as the development of numerical schemes based on the dual variable method and the Lions transversal variable method [1].

## Funding statement

D.G. was partially supported by KAUST SRI, Uncertainty Quantification Center in Computational Science and Engineering and CAMGSD-LARSys (FCT-Portugal), PTDC/MAT-CAL/0749/ 2012 (FCT-Portugal). R.M.V. was partially supported by CNPq-Brazil through a PhD scholarship ‘Programme Science without Borders and KAUST’ Saudi Arabia. M.-T.W. acknowledges support from the Austrian Academy of Sciences ÖAW via the New Frontiers Project NST-001.

## Footnotes

One contribution of 13 to a Theme Issue ‘Partial differential equation models in the socio-economic sciences’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.