*population* of actors, each of which has a level of wealth (a single number) that changes over time. On each time step two agents (chosen by an *interaction rule*) interact with each other and exchange wealth (according to a *transaction rule*). The idea is to understand the evolution of the population's wealth over time. My hazy memory is that this idea came from a class by Prof. Sven Anderson at Bard (any errors or misconceptions here are due to my (Peter Norvig) misunderstanding of his idea). Why this is interesting: (1) an example of using simulation to model the world. (2) Many students will have preconceptions about how economies work that will be challenged by the results shown here.

`sample`

function, which samples *N* elements form a distribution and then normalizes them to have a given mean. By default we will have *N*=5000 actors and an initial mean wealth of 100 simoleons.

In [299]:

```
import random
import matplotlib
import matplotlib.pyplot as plt
N = 5000 # Default size of population
mu = 100. # Default mean of population's wealth
def sample(distribution, N=N, mu=mu):
"Sample from the distribution N times, then normalize results to have mean mu."
return normalize([distribution() for _ in range(N)], mu * N)
def constant(mu=mu): return mu
def uniform(mu=mu, width=mu): return random.uniform(mu-width/2, mu+width/2)
def gauss(mu=mu, sigma=mu/3): return random.gauss(mu, sigma)
def beta(alpha=2, beta=3): return random.betavariate(alpha, beta)
def pareto(alpha=4): return random.paretovariate(alpha)
def normalize(numbers, total):
"Scale the numbers so that they add up to total."
factor = total / float(sum(numbers))
return [x * factor for x in numbers]
```

In [360]:

```
def random_split(X, Y):
"Take all the money in the pot and divide it randomly between X and Y."
pot = X + Y
m = random.uniform(0, pot)
return m, pot - m
def winner_take_most(X, Y, most=3/4.):
"Give most of the money in the pot to one of the parties."
pot = X + Y
m = random.choice((most * pot, (1 - most) * pot))
return m, pot - m
def winner_take_all(X, Y):
"Give all the money in the pot to one of the actors."
return winner_take_most(X, Y, 1.0)
def redistribute(X, Y):
"Give 55% of the pot to the winner; 45% to the loser."
return winner_take_most(X, Y, 0.55)
def split_half_min(X, Y):
"""The poorer actor only wants to risk half his wealth;
the other actor matches this; then we randomly split the pot."""
pot = min(X, Y)
m = random.uniform(0, pot)
return X - pot/2. + m, Y + pot/2. - m
```

`anyone`

samples two members of the population uniformly and independently, but there are other possible rules, like `nearby(pop, k)`

, which choses one member uniformly and then chooses a second within `k`

index elements away, to simulate interactions within a local neighborhood.

In [356]:

```
def anyone(pop): return random.sample(range(len(pop)), 2)
def nearby(pop, k=5):
i = random.randrange(len(pop))
j = i + random.choice((1, -1)) * random.randint(1, k)
return i, (j % len(pop))
def nearby1(pop): return nearby(pop, 1)
```

Now let's describe the code to run the simulation and summarize/plot the results. The function `simulate`

does the work; it runs the interaction function to find two actors, then calls the transaction function to figure out how to split their wealth, and repeats this T times. The only other thing it does is record results. Every so-many steps, it records some summary statistics of the population (by default, this will be every 25 steps).

What information do we record to summarize the population? Out of the N=5000 (by default) actors, we will record the wealth of exactly nine of them: the ones, in sorted-by-wealth order that occupy the 1% spot (that is, if N=5000, this would be the 50th wealthiest actor), then the 10%, 25% 1/3, and median; and then likewise from the bottom the 1%, 10%, 25% and 1/3.

(Note that we record the *median*, which changes over time; the *mean* is defined to be 100 when we start, and since all transactions conserve wealth, the mean will always be 100.)

What do we do with these results, once we have recorded them? First we print them in a table for the first time step, the last, and the middle. Then we plot them as nine lines in a plot where the y-axis is wealth and the x-axis is time (note that when the x-axis goes from 0 to 1000, and we have `record_every=25`

, that means we have actually done 25,000 transactions, not 1000).

In [368]:

```
def simulate(population, transaction_fn, interaction_fn, T, percentiles, record_every):
"Run simulation for T steps; collect percentiles every 'record_every' time steps."
results = []
for t in range(T):
i, j = interaction_fn(population)
population[i], population[j] = transaction_fn(population[i], population[j])
if t % record_every == 0:
results.append(record_percentiles(population, percentiles))
return results
def report(distribution=gauss, transaction_fn=random_split, interaction_fn=anyone, N=N, mu=mu, T=5*N,
percentiles=(1, 10, 25, 33.3, 50, -33.3, -25, -10, -1), record_every=25):
"Print and plot the results of the simulation running T steps."
# Run simulation
population = sample(distribution, N, mu)
results = simulate(population, transaction_fn, interaction_fn, T, percentiles, record_every)
# Print summary
print('Simulation: {} * {}(mu={}) for T={} steps with {} doing {}:\n'.format(
N, name(distribution), mu, T, name(interaction_fn), name(transaction_fn)))
fmt = '{:6}' + '{:10.2f} ' * len(percentiles)
print(('{:6}' + '{:>10} ' * len(percentiles)).format('', *map(percentile_name, percentiles)))
for (label, nums) in [('start', results[0]), ('mid', results[len(results)//2]), ('final', results[-1])]:
print fmt.format(label, *nums)
# Plot results
for line in zip(*results):
plt.plot(line)
plt.show()
def record_percentiles(population, percentiles):
"Pick out the percentiles from population."
population = sorted(population, reverse=True)
N = len(population)
return [population[int(p*N/100.)] for p in percentiles]
def percentile_name(p):
return ('median' if p == 50 else
'{} {}%'.format(('top' if p > 0 else 'bot'), abs(p)))
def name(obj):
return getattr(obj, '__name__', str(obj))
```

Finally, let's run a simulation!

In [369]:

```
report(gauss, random_split)
```

In [376]:

```
report(uniform, random_split)
report(beta, random_split)
report(pareto, random_split)
report(constant, random_split)
```

`random_split`

rule then in the end, wealth accumulates to the top third at the expense of the bottom two-thirds, regardless of starting population.

`random_split`

rule produces inequality: the actor at the bottom quarter of the population has only about a third of the mean wealth, and the actor at the top 1% spot has 4.5 times the mean.
Suppose we want a society with more income equality. We could use the `split_half_min`

rule, in which each transaction has a throttle in that the poorer party only risks half of their remaining wealth. Or we could use the `redistribute`

rule, in which the loser of a transaction still gets 45% of the total (meaning the loser will actually gain in many transactions). Let's see what effects these rules have. In analyzing these plots, note that they have different Y-axes.

In [371]:

```
report(gauss, random_split)
report(gauss, redistribute)
report(gauss, split_half_min)
```

`redistribute`

rule is very effective in reducing income inequality: the lines of the plot all converge towards the mean of 100 instead of diverging. With the `split_half_min`

rule, inequality increases at a rate about half as fast as `random_split`

. However, the `split_half_min`

plot looks like it hasn't converged yet (whereas all the other plots reach convergence at about the 500 mark). Let's try running `split_half_min`

10 times longer:

In [372]:

```
report(gauss, split_half_min, T=50*N)
```

It looks like `split_half_min`

*still* hasn't converged, and is continuing to (slowly) drive wealth to the top 10%.

Now let's shift gears: suppose that we don't care about decreasing income inequality; instead we want to increase opportunity for some actors to become wealthier. We can try the `winner_take_most`

or `winner_take_all`

rules (compared to the baseline `random_split`

):

In [373]:

```
report(gauss, random_split)
report(gauss, winner_take_most)
report(gauss, winner_take_all)
```

We see that the `winner_take_most`

rule, in which the winner of a transaction takes 3/4 of the pot, does not increase the opportunity for wealth as much as `random_split`

, but that `winner_take_all`

is very effective at concentrating almost all the wealth in the hands of the top 10%, and makes the top 1% 4 times as wealthy as `random_split`

.

That suggests we look at where the breaking point is. Let's consider several different amounts for what winner takes:

In [375]:

```
def winner_take_80(X, Y): return winner_take_most(X, Y, 0.80)
def winner_take_90(X, Y): return winner_take_most(X, Y, 0.90)
def winner_take_95(X, Y): return winner_take_most(X, Y, 0.95)
report(gauss, winner_take_80)
report(gauss, winner_take_90)
report(gauss, winner_take_95)
```

`random_split`

, and that winner takes 95% is similar to winner takes all for the top 10%, but is much kinder to the bottom 75%.

*local*; that you can only do business with your close neighbors. Will that make income more equitable, because there will be no large, global conglomorates? Let's see:

In [374]:

```
report(gauss, random_split, anyone)
report(gauss, random_split, nearby)
report(gauss, random_split, nearby1)
```

`nearby`

rule, which limits transactions to your 5 closest neighbors in either direction (out of 5000 total actors), has a negligible effect on the outcome. I found that fairly surprising. But the `nearby1`

rule, which lets you do business only with your immediate left or right neighbor does have a slight effect towards income equality. The bottom quarter still do poorly, but the top 1% only gets to about 85% of what they get under unconstrained trade.

In [ ]:

```
```