Wealth Inequality in an Exchange Economy

def gini(wealth):
    """Gini coefficient of a wealth distribution.

    Uses the sorted-array formula:
        G = (2 * sum(i * w_(i)) - (n+1) * sum(w)) / (n * sum(w))
    where indices i are 1-based and w_(i) is the i-th smallest wealth.
    Returns 0 for perfect equality and (n-1)/n when one agent holds everything.
    """
    n = len(wealth)
    sorted_w = np.sort(wealth)
    indices = np.arange(1, n + 1)
    total = np.sum(sorted_w)
    return (2.0 * np.dot(indices, sorted_w) - (n + 1) * total) / (n * total)

def simulate_exchange(n_agents=N_AGENTS, n_steps=N_STEPS, seed=SEED):
    """Run the random-exchange wealth model.

    Agents start with equal wealth of 1.0 each.  At every step two distinct
    agents are chosen uniformly at random; the poorer transfers a uniformly
    random fraction of their own wealth to the richer.  Total wealth is
    conserved exactly at n_agents throughout.

    Returns (final_wealth, gini_history) where gini_history[t] is the Gini
    coefficient after t exchanges (gini_history[0] is the initial value 0.0).
    """
    rng = np.random.default_rng(seed)
    wealth = np.ones(n_agents, dtype=float)
    gini_history = np.empty(n_steps + 1)
    gini_history[0] = gini(wealth)
    for step in range(n_steps):
        i, j = rng.choice(n_agents, size=2, replace=False)
        if wealth[i] > wealth[j]:
            i, j = j, i  # i is now the poorer agent
        fraction = rng.uniform(0.0, 1.0)
        transfer = fraction * wealth[i]
        wealth[i] -= transfer
        wealth[j] += transfer
        gini_history[step + 1] = gini(wealth)
    return wealth, gini_history

def lorenz_curve(wealth):
    """Return (population_fractions, wealth_fractions) arrays for the Lorenz curve.

    Both arrays start at (0, 0) and end at (1, 1).  The area between the
    curve and the 45-degree equality line equals G/2, where G is the Gini
    coefficient.
    """
    n = len(wealth)
    sorted_w = np.sort(wealth)
    cum_w = np.cumsum(sorted_w)
    pop_fracs = np.arange(0, n + 1) / n
    wealth_fracs = np.concatenate([[0.0], cum_w / cum_w[-1]])
    return pop_fracs, wealth_fracs

def plot_lorenz_curve(wealth, filename):
    """Save a Lorenz curve with the line of equality overlaid."""
    pop_fracs, wealth_fracs = lorenz_curve(wealth)
    df_lorenz = pl.DataFrame({"population": pop_fracs, "wealth": wealth_fracs})
    df_equal = pl.DataFrame({"population": [0.0, 1.0], "wealth": [0.0, 1.0]})
    lorenz_line = (
        alt.Chart(df_lorenz)
        .mark_line(color="steelblue", strokeWidth=2)
        .encode(
            x=alt.X("population:Q", title="Cumulative share of population"),
            y=alt.Y("wealth:Q", title="Cumulative share of wealth"),
        )
    )
    equal_line = (
        alt.Chart(df_equal)
        .mark_line(color="gray", strokeDash=[4, 4], strokeWidth=1)
        .encode(x="population:Q", y="wealth:Q")
    )
    chart = alt.layer(equal_line, lorenz_line).properties(
        width=350, height=350, title="Lorenz curve after 2000 exchanges"
    )
    chart.save(filename)

def plot_gini_trajectory(gini_history, filename):
    """Save a line chart of the Gini coefficient at each exchange step."""
    df = pl.DataFrame(
        {
            "step": np.arange(len(gini_history)),
            "gini": gini_history,
        }
    )
    chart = (
        alt.Chart(df)
        .mark_line(color="steelblue", strokeWidth=2)
        .encode(
            x=alt.X("step:Q", title="Exchange step"),
            y=alt.Y(
                "gini:Q", scale=alt.Scale(domain=[0.0, 1.0]), title="Gini coefficient"
            ),
        )
        .properties(width=450, height=300, title="Gini coefficient over time")
    )
    chart.save(filename)

import numpy as np
import pytest
from wealth import gini, simulate_exchange, N_AGENTS


def test_gini_equal_distribution():
    # With identical wealth for every agent the Gini must be exactly 0.
    # Analytically: (2 * n(n+1)/2 - (n+1)*n) / n^2 = 0 for any n.
    assert gini(np.ones(50)) == pytest.approx(0.0, abs=1e-12)


def test_gini_perfect_inequality():
    # One agent holds all wealth; Gini = (n-1)/n.
    # With n=10: G = 9/10 = 0.9 exactly.
    n = 10
    w = np.zeros(n)
    w[-1] = float(n)
    assert gini(w) == pytest.approx((n - 1) / n, rel=1e-9)


def test_gini_bounded_during_simulation():
    # Wealth is always non-negative and total wealth is conserved, so the
    # Gini coefficient must stay in [0, 1] at every step.
    _, gini_history = simulate_exchange()
    assert np.all(gini_history >= 0.0)
    assert np.all(gini_history <= 1.0)


def test_total_wealth_conserved():
    # Every exchange transfers wealth without creation or destruction, so
    # the sum must equal n_agents (each agent starts with 1.0) to within
    # floating-point rounding.
    final_wealth, _ = simulate_exchange()
    assert np.sum(final_wealth) == pytest.approx(N_AGENTS, rel=1e-9)


def test_gini_rises_from_zero():
    # The model starts from perfect equality (Gini = 0) and should
    # produce substantial inequality after 2000 exchanges.
    # The theoretical steady-state Gini for exponential wealth distributions
    # is 0.5; after 2000 steps the simulated value reliably exceeds 0.3.
    _, gini_history = simulate_exchange()
    assert gini_history[0] == pytest.approx(0.0, abs=1e-12)
    assert gini_history[-1] > 0.3