This is a simulation of an economic marketplace in which there is a 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.
First things first: what should our initial population look like? We will provide several distribution functions (constant, uniform, Gaussian, etc.) and a 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 a transaction, two actors come together; they have existing wealth levels X and Y. For now we will only consider transactions that conserve wealth, so our transaction rules will decide how to split up the pot of X+Y total wealth.
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
How do you decide which parties interact with each other? The rule 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)
How do we interpret this? Well, we can see the mass of wealth spreading out: the rish get richer and the poor get poorer. We know the rich get richer because the blue and green lines (top 10% and top 1%) are going up: the actor in the 1% position (the guy with the least money out of the 1%, or to put it another way, the most money out of the 99%) starts with 177.13 and ends up with 447.98 (note this is not necessarily the same guy, just the guy who ends up in that position). The guy at the 10% spot also gets richer, going from 141.87 to 228.06. The 25% and 33% marks stay roughly flat, but everyone else gets poorer! The median actor loses 30% of his wealth, and the bottom 1% actor loses almost 95% of his wealth.
Now let's see if the starting population makes any difference. My vague impression is that we're dealing with ergodic Markov chains and it doesn't much matter what state you start in. But let's see:
In [376]:
report(uniform, random_split)
report(beta, random_split)
report(pareto, random_split)
report(constant, random_split)
It looks like we can confirm that the starting population doesn't matter much—if we are using the 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.
Now let's see what happens when we vary the transaction rule. The 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)
We see that the 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)
We see that winner takes 80% produces results similar to 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%.
Suppose that transactions are constrained to be 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)
We see that the 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 [ ]: