Level 6 · Economics: Dynamics and General Equilibrium

Until now each market was studied on its own: one price, one supply curve, one demand curve, everything else held still. The micro half of this level lets every market move at once. In an exchange economy people's incomes are the value of what they own, so a change in one price changes everybody's income and every demand. A general equilibrium is a list of prices at which every market clears together. You will find it for two goods and for many, watch an imaginary auctioneer grope towards it, ask whether the result could be improved for anyone without hurting someone else (the first welfare theorem says no), and follow orders as they ripple through industries that supply each other (the input-output model). Most of it comes down to small linear systems, solved by hand.

The macro half adds time. People choose today knowing that today's choice changes what they can do tomorrow: how much to spend and how much to save, when to eat the stores, which job offer to take. The tools are the Euler equation ("moving one coin from today to tomorrow cannot help") and the Bellman equation, which values a situation by the best choice in it plus the value of where that choice leads. Then come markets whose dynamics are set by simple rules instead of foresight: farmers who plant on last year's price, customers who follow a queue, traders who meet at random.

The paths meet in exact places. Backward induction and memoised recursion are the dynamic programming of the Algorithms path (Level 5). Gauss-Jordan elimination, determinants and the matrix-vector product are the matrices of the Quantum Computing path (Level 1), and the linear systems of the Signals & Control path. Value function iteration is a fixed-point iteration that shrinks the error by the same factor every step, like the Riccati iteration of the Control path's LQR (which is itself a Bellman equation), and the cobweb model is a first-order difference equation whose gap is multiplied by a negative factor each step. Agent-based markets are seeded simulations, exactly as on the Statistics path.

Exchange: income is what you own, at today's prices

In an exchange economy nobody earns a wage. Each household i owns a stock of goods, its endowment: e_{ij} units of each good j. It comes to the market to trade. At prices p_j selling everything would raise

m_i = \sum_j p_j\, e_{ij}

so that is the household's income: it can afford any bundle worth no more than its endowment. With Level 2's Cobb-Douglas tastes the household spends the fixed share s_{ij} of its income on good j, so it demands

x_{ij} = \frac{s_{ij}\, m_i}{p_j}

When tastes are given as positive weights w_{ij}, the shares are the weights divided by their total, s_{ij} = w_{ij} / \sum_k w_{ik}, as in Level 2. Two things are new. Income now depends on the prices: a higher price of grain makes a grain farmer richer. And if every price is multiplied by the same factor, every income is multiplied by it too, so no demand changes. Only relative prices matter.

What a market cares about is the excess demand for each good: what all households want minus what they own between them,

z_j = \sum_i x_{ij} - \sum_i e_{ij}

Positive means the good is short at these prices; negative means some is left over. A potter and a farmer trade bowls, cloth and grain:

weights = [[1, 2, 1], [3, 1, 4]]          # tastes as weights, one row per household
owns = [[12, 0, 10], [0, 6, 20]]          # bowls, cloth, grain
total = [sum(e[j] for e in owns) for j in range(3)]

def demands(prices):
    out = []
    for w, e in zip(weights, owns):
        m = sum(p * q for p, q in zip(prices, e))                  # income: the endowment's value
        out.append([w[j] / sum(w) * m / prices[j] for j in range(3)])
    return out

for prices in ([2, 5, 1], [4, 10, 2], [2, 3, 1]):
    x = demands(prices)
    z = [sum(row[j] for row in x) - total[j] for j in range(3)]
    print(prices, " demands", [[round(v, 3) for v in row] for row in x])
    print("   excess", [round(v, 3) + 0.0 for v in z], " value", round(sum(p * zj for p, zj in zip(prices, z)), 9) + 0.0)

[2, 5, 1]  demands [[4.25, 3.4, 8.5], [9.375, 1.25, 25.0]]
   excess [1.625, -1.35, 3.5]  value 0.0
[4, 10, 2]  demands [[4.25, 3.4, 8.5], [9.375, 1.25, 25.0]]
   excess [1.625, -1.35, 3.5]  value 0.0
[2, 3, 1]  demands [[4.25, 5.667, 8.5], [7.125, 1.583, 19.0]]
   excess [-0.625, 1.25, -2.5]  value 0.0

Doubling every price changed nothing. None of these prices clears the markets: at (2, 5, 1) bowls and grain are short and cloth is left over, at (2, 3, 1) the reverse. The last line of each pair is Walras' law: the value of the excess demands, \sum_j p_j z_j, is zero at every price vector, cleared markets or not. The reason is a budget: each household spends exactly its income, so the value of what it buys equals the value of what it owns, \sum_j p_j x_{ij} = \sum_j p_j e_{ij}. Add this over the households and move everything to one side. In floats the sum comes out as a tiny number such as 1e-15, which is rounding; a value far from zero means a demand is wrong, so Walras' law makes a good test of your own code.

Common mistakes: treating the endowment as what the household consumes (it trades away from it). Using raw weights as shares. Keeping income fixed when prices change. Forgetting to subtract the endowments, which gives total demand instead of excess demand. Expecting \sum_j p_j z_j to be exactly 0.0 in floats.

Practice: Wool for Apples, The Harbour Ledger Always Balances

One price is free: the numéraire and the clearing price

Since only relative prices matter, fix one of them. The good whose price is set to 1 is the numéraire: every other price is then "how many units of the numéraire one unit of this good trades for". With good 1 as the numéraire, p_1 = 1, a two-good economy has one unknown left, p_2.

A competitive equilibrium is a price at which every market clears: z_j = 0 for all j. With two goods, clearing one market is enough. Walras' law says p_1 z_1 + p_2 z_2 = 0, so if z_2 = 0 then p_1 z_1 = 0, and z_1 = 0 because p_1 = 1. With m goods, clearing any m - 1 markets clears the last one.

Three households trade salt (good 1) and honey (good 2); shares[i] is household i's share of income spent on salt. The excess demand for honey falls as honey gets dearer, so Level 1's bisection finds the price where it crosses zero:

shares = [0.7, 0.2, 0.5]
owns = [(10, 2), (0, 15), (4, 4)]               # (salt, honey)
E1, E2 = sum(e1 for e1, e2 in owns), sum(e2 for e1, e2 in owns)

def excess(p2):                                 # (z1, z2) at prices (1, p2)
    z1, z2 = -E1, -E2
    for a, (e1, e2) in zip(shares, owns):
        m = e1 + p2 * e2
        z1 += a * m
        z2 += (1 - a) * m / p2
    return z1, z2

for p2 in (0.25, 0.5, 1, 2):
    print("p2 =", p2, " excess demands", [round(v, 4) for v in excess(p2)])

lo, hi = 1e-6, 1e6                              # honey short at lo, left over at hi
while hi - lo > 1e-12 * hi:
    mid = (lo + hi) / 2
    if excess(mid)[1] > 0:
        lo = mid
    else:
        hi = mid
print("clearing price", round(lo, 6), " excess there", [round(v, 9) + 0.0 for v in excess(lo)])

p2 = 0.25  excess demands [-3.4, 13.6]
p2 = 0.5  excess demands [-1.8, 3.6]
p2 = 1  excess demands [1.4, -1.4]
p2 = 2  excess demands [7.8, -3.9]
clearing price 0.78125  excess there [0.0, 0.0]

At 0.25 honey is short and salt is left over; at 2 the reverse. The bisection stops at a relative width, since prices can be of any size. Both markets clear at the price it finds, although only the honey market was used. Press Run to see the honey market's excess demand cross zero once:

shares, owns = [0.7, 0.2, 0.5], [(10, 2), (0, 15), (4, 4)]
ps = [k / 100 for k in range(20, 301)]
z2 = [sum((1 - a) * (e1 + p * e2) / p for a, (e1, e2) in zip(shares, owns)) - 21 for p in ps]
plot(ps, z2, label="excess demand for honey", xlabel="price of honey (salt = 1)")
plot(ps, [0 for p in ps], label="cleared")

Bisection is not the only way. Multiply the honey market's clearing equation by p_2, and income, which is linear in p_2, makes the whole equation linear: one unknown, solved by moving terms across. Doing that for m goods at once is the subject of a later section.

Common mistakes: forgetting that income moves with p_2 (the honey a household owns is worth more when honey is dear). Leaving p_1 free, so that any multiple of the answer also "works". Solving for p_1 / p_2 when p_2 / p_1 is asked for. Stopping a bisection on an absolute width when the price may be tiny or huge.

Practice: What a Kilo of Apples Fetches

Tâtonnement: an auctioneer feels for the prices

Nobody in a real market knows the equilibrium prices in advance. Léon Walras pictured an auctioneer who calls out prices, hears how much of each good people would want, raises the price of what is short and lowers the price of what is left over, and lets trade happen only once every market clears. This groping is called tâtonnement: Level 1's price-adjustment loop, run for every market at once, with each market's demand moving whenever any price moves.

A common rule moves each price in proportion to its relative excess demand, z_j / E_j with E_j the total that exists, so that the same step means the same thing for a good the economy has 3 of and one it has 3000 of:

p_j \leftarrow p_j \left(1 + \text{step} \cdot \frac{z_j}{E_j}\right)

Nobody can want less than nothing, so z_j \ge -E_j, and with a step below 1 every price stays positive. Problems state their exact loop, and the order of its parts matters: measure the excess demands, test the tolerance, give up if the budget of rounds is spent, adjust every price from the same measurement, and, if the numéraire moved too, divide every price by its new value. For the salt and honey market, holding salt at 1 and moving only honey:

shares, owns = [0.7, 0.2, 0.5], [(10, 2), (0, 15), (4, 4)]
E2 = sum(e2 for e1, e2 in owns)

def z2(p2):
    return sum((1 - a) * (e1 + p2 * e2) / p2 for a, (e1, e2) in zip(shares, owns)) - E2

for step in (0.5, 2, 4):
    p2, path = 2.0, []
    for rnd in range(8):
        path.append(round(p2, 4))
        p2 = p2 * (1 + step * z2(p2) / E2)
    print("step", step, path)

step 0.5 [2.0, 1.8143, 1.6569, 1.5234, 1.4103, 1.3145, 1.2332, 1.1644]
step 2 [2.0, 1.2571, 0.9671, 0.8538, 0.8096, 0.7923, 0.7856, 0.7829]
step 4 [2.0, 0.5143, 0.8397, 0.7684, 0.7841, 0.7806, 0.7814, 0.7812]

The clearing price is 0.78125 (the previous section found it). With step 0.5 the calls crawl towards it from one side; with step 2 they get there much faster; with step 4 each correction overshoots, and the calls zigzag around the price before they settle. The auctioneer is a feedback loop, and as on the Signals & Control path, a gain that is too high makes a loop swing; higher still, it can diverge, or push a price below zero. For some economies (with tastes other than Cobb-Douglas) tâtonnement never settles at any step size, which is one reason to also know the exact methods later in this lesson.

Common mistakes: testing the tolerance after adjusting instead of before. Updating p[0] first and then using the new value for the other goods. Renormalising by the old numéraire price. Using the absolute excess demand z_j where the rule says z_j / E_j. Returning "gave up" one round early or late.

Practice: The Auctioneer Calls Out Prices

The Edgeworth box: who should get more?

Now forget prices and ask what a good split looks like. Two people, A and B, share fixed totals X of good 1 and Y of good 2. Every split is a point in a rectangle of width X and height Y, the Edgeworth box: if A has (x, y), B has the rest, (X - x, Y - y).

A split is Pareto efficient if no other split makes one person better off without making the other worse off. A split that does is a Pareto improvement, and it dominates the first: nobody worse off, somebody strictly better off.

The key measure is the marginal rate of substitution (MRS): how much of good 2 a person would give up for a tiny extra amount of good 1, the ratio of the two marginal utilities. Level 2 found the marginal utilities of u = x^a y^{1-a}, namely a\,u/x and (1-a)\,u/y, so

\text{MRS} = \frac{a}{1 - a} \cdot \frac{y}{x}

Inside the box (both have some of each good), a split is efficient exactly when the two MRS are equal. If A's MRS is larger, A values good 1 more than B does, in terms of good 2. Then A buys a little good 1 from B and pays in good 2 at any rate between the two MRS: A gets good 1 for less than it is worth to A, and B is paid more than it is worth to B. Both gain. Bread (good 1) and cheese shared by two walkers:

a, b, X, Y = 0.3, 0.6, 12, 20

def u(s, x, y):
    return x ** s * y ** (1 - s)

def mrs(s, x, y):
    return s / (1 - s) * y / x

x, y = 4, 10                                    # A's bundle; B has (8, 10)
h = 1e-6
steps = ((u(a, x + h, y) - u(a, x, y)) / h) / ((u(a, x, y + h) - u(a, x, y)) / h)
print("A's MRS", round(mrs(a, x, y), 4), " from tiny steps", round(steps, 4), "  B's MRS", round(mrs(b, X - x, Y - y), 4))

dx, rate = 0.1, 1.5                             # B buys 0.1 bread and pays 1.5 cheese per bread
before = (u(a, x, y), u(b, X - x, Y - y))
after = (u(a, x - dx, y + rate * dx), u(b, X - x + dx, Y - y - rate * dx))
print("before", [round(v, 4) for v in before], " after", [round(v, 4) for v in after])

A's MRS 1.0714  from tiny steps 1.0714   B's MRS 1.875
before [7.5966, 8.7469]  after [7.6181, 8.7592]

B's MRS is the larger, so B should get more bread, and a trade at 1.5 cheese per bread (between 1.07 and 1.875) makes both walkers better off: the split (4, 10) was not efficient.

The contract curve is the set of all efficient splits. For each amount of bread x that A might hold, find the cheese y at which the two MRS are equal. A's MRS rises with y and B's falls, so their difference crosses zero once and bisection finds it (the equation can also be solved by algebra: it is linear in y once the fractions are multiplied out).

a, b, X, Y = 0.3, 0.6, 12, 20

def gap(x, y):                                  # A's MRS minus B's
    return a / (1 - a) * y / x - b / (1 - b) * (Y - y) / (X - x)

curve = []
for x in (2, 4, 6, 8, 10):
    lo, hi = 0.0, Y
    for _ in range(100):
        mid = (lo + hi) / 2
        if gap(x, mid) < 0:
            lo = mid
        else:
            hi = mid
    curve.append((x, round(lo, 4)))
print(curve)

[(2, 8.2353), (4, 12.7273), (6, 15.5556), (8, 17.5), (10, 18.9189)]

At x = 4 the efficient amount of cheese for A is far above 10: from (4, 10), the trade above moves towards the curve. The curve runs from one corner of the box to the other. Press Run to draw it:

a, b, X, Y = 0.3, 0.6, 12, 20
ka, kb = a / (1 - a), b / (1 - b)
xs = [X * k / 100 for k in range(1, 100)]
ys = []
for x in xs:
    lo, hi = 0.0, Y
    for _ in range(60):
        mid = (lo + hi) / 2
        lo, hi = (mid, hi) if ka * mid / x < kb * (Y - mid) / (X - x) else (lo, mid)
    ys.append(lo)
plot(xs, ys, label="contract curve", xlabel="A's bread", ylabel="A's cheese", title="Edgeworth box")

Two floats computed in different ways rarely agree to the last digit, so "equal MRS" needs a tolerance, usually a relative one, |\text{MRS}_A - \text{MRS}_B| \le 10^{-9} \max(\text{MRS}_A, \text{MRS}_B).

Common mistakes: putting A's bundle into B's MRS (B's MRS uses (X - x, Y - y)). Swapping the labels: the person with the larger MRS should get more of good 1. Comparing MRS with ==. Evaluating the MRS at the corners of the box, where a zero divides.

Practice: Splitting the Trail Supplies

The first welfare theorem

Put prices back. In a competitive equilibrium each household picks its best affordable bundle, and the best Cobb-Douglas bundle sets the MRS equal to the price ratio p_1 / p_2 (Level 2's rule of equal marginal utility per coin, rewritten). Every household faces the same prices, so every household ends with the same MRS: the equilibrium is on the contract curve. That is the first welfare theorem: a competitive equilibrium allocation is Pareto efficient. It holds under stated assumptions: everyone takes prices as given, there are no externalities (Level 5), there is a market for every good, and more is always better.

Check it numerically for the walkers, now with endowments (10, 4) for A and (2, 16) for B, by searching a grid of splits for one that both prefer. Comparing utilities also needs a tolerance, here |u - v| \le 10^{-9} \max(1, |u|, |v|), so write one comparison function and build every verdict from it:

a, b = 0.3, 0.6
eA, eB = (10, 4), (2, 16)
X, Y = eA[0] + eB[0], eA[1] + eB[1]

def z2(p2):                                         # excess demand for cheese, bread = 1
    mA, mB = eA[0] + p2 * eA[1], eB[0] + p2 * eB[1]
    return ((1 - a) * mA + (1 - b) * mB) / p2 - Y

lo, hi = 1e-6, 1e6
while hi - lo > 1e-13 * hi:
    mid = (lo + hi) / 2
    lo, hi = (mid, hi) if z2(mid) > 0 else (lo, mid)
p2 = lo
mA = eA[0] + p2 * eA[1]
xA, yA = a * mA, (1 - a) * mA / p2
uA, uB = xA ** a * yA ** (1 - a), (X - xA) ** b * (Y - yA) ** (1 - b)
print("price", round(p2, 6), " A gets", (round(xA, 4), round(yA, 4)), " MRS", round(a / (1 - a) * yA / xA, 6), round(b / (1 - b) * (Y - yA) / (X - xA), 6), " 1/p2", round(1 / p2, 6))

def cmp(u, v):                                      # 1: u better, -1: v better, 0: equal
    if abs(u - v) <= 1e-9 * max(1.0, abs(u), abs(v)):
        return 0
    return 1 if u > v else -1

count = {}
for i in range(0, 25):
    for j in range(0, 41):
        x, y = i / 2, j / 2
        key = (cmp(x ** a * y ** (1 - a), uA), cmp((X - x) ** b * (Y - y) ** (1 - b), uB))
        count[key] = count.get(key, 0) + 1
print(sorted(count.items()))

price 0.722222  A gets (3.8667, 12.4923)  MRS 1.384615 1.384615  1/p2 1.384615
[((-1, -1), 239), ((-1, 1), 370), ((1, -1), 416)]

The two MRS agree, and both equal the price ratio. In the count, a key (1, 1) would be a split both walkers strictly prefer; it never appears, and neither does a split that is better for one and equal for the other. Splits better for one walker exist (keys (1, -1) and (-1, 1)), but each costs the other. The utility of a bundle with none of one good is 0 (0 ** 0.3 == 0.0), so the edges of the box need no special case.

The theorem says nothing about fairness: an equilibrium from very unequal endowments is efficient too. Which efficient point is reached depends on who owned what at the start.

Aside (optional): what the welfare theorems are used for. Some economists read the first welfare theorem as a strong argument for markets: under its assumptions decentralised trade wastes nothing, with no planner needing to know anyone's tastes. Others stress how demanding the assumptions are (no market power, no externalities, complete markets, full information) and read it as a checklist of the ways real markets can fail, each a reason to look closely at a market. Large computable general equilibrium models, with dozens of industries and households, are used by governments and international bodies to estimate the effects of trade agreements and taxes. Supporters value that they keep every market and budget consistent at once; critics point out that results depend heavily on assumed elasticities and on the model's structure, which data pin down only loosely. Problems here state their model and ask what follows in it.

Common mistakes: testing "one is better" before "both are equal" (two equal comparisons pass both dominance tests). Comparing utilities with ==. Giving B the proposal itself instead of the rest of the box. Reading "efficient" as "fair".

Practice: The Market's Split Is Hard to Beat

Industries that supply each other: input-output

Exchange economies have no production. Real industries buy from each other: the farms buy tools, the tool makers buy transport, the carriers buy grain for their horses. Wassily Leontief's input-output model writes this as a matrix. A_{ij} is how many dollars of industry i's output are used up to make one dollar of industry j's output, so column j is industry j's recipe, and row i lists what industry i supplies to each buyer. On top of that, households, visitors and the government want d_i dollars of each industry's output: the final demand.

Total output x must cover the final demand and everything the industries use up:

x = A\,x + d, \qquad (A x)_i = \sum_j A_{ij}\, x_j

Think of it as orders spreading. To deliver d, industries need inputs A d; to make those inputs they need A(Ad) = A^2 d; and so on. Total output is the sum of the layers

x = d + A d + A^2 d + A^3 d + \cdots

exactly like Level 2's rounds of re-spending behind the multiplier. Each layer is smaller than the one before when every column of A sums to less than 1 (every industry adds value of its own): the total of a layer is \sum_j (\text{column sum}_j) \times (\text{previous layer})_j, at most the largest column sum times the previous total. Three industries, grain, tools and transport:

A = [[0.10, 0.20, 0.05],          # A[i][j]: dollars of industry i's output per dollar of industry j's
     [0.15, 0.10, 0.30],
     [0.20, 0.25, 0.10]]
d = [50, 30, 20]
n = len(d)
print("column sums", [round(sum(A[i][j] for i in range(n)), 2) for j in range(n)])

def times(A, v):                  # the matrix-vector product A v
    return [sum(A[i][j] * v[j] for j in range(len(v))) for i in range(len(A))]

layer, total = d, list(d)
for k in range(1, 6):
    new = times(A, layer)
    print("layer", k, [round(v, 3) for v in new], " size", round(sum(new), 3), " ratio", round(sum(new) / sum(layer), 3))
    layer = new
    total = [t + v for t, v in zip(total, layer)]
print("layers 0 to 5", [round(v, 3) for v in total])
for k in range(6, 200):
    layer = times(A, layer)
    total = [t + v for t, v in zip(total, layer)]
print("layers 0 to 199", [round(v, 4) for v in total], " check x - Ax - d:", [round(v, 9) + 0.0 for v in [t - a - b for t, a, b in zip(total, times(A, total), d)]])

column sums [0.45, 0.55, 0.45]
layer 1 [12.0, 16.5, 19.5]  size 48.0  ratio 0.48
layer 2 [5.475, 9.3, 8.475]  size 23.25  ratio 0.484
layer 3 [2.831, 4.294, 4.268]  size 11.393  ratio 0.49
layer 4 [1.355, 2.134, 2.066]  size 5.556  ratio 0.488
layer 5 [0.666, 1.037, 1.011]  size 2.714  ratio 0.488
layers 0 to 5 [72.327, 63.265, 55.32]
layers 0 to 199 [72.9614, 64.2551, 56.2845]  check x - Ax - d: [0.0, 0.0, 0.0]

Each layer is about half the size of the last, never more than 0.55 of it (the largest column sum), and a final demand of 100 dollars needs about 193.5 dollars of output in all. The layers are also what the update x \leftarrow d + A x produces, started from x = d: after k updates x is the sum of layers 0 to k. That iteration needs no matrix algebra at all, only a stopping rule, which problems state. Build each new vector entirely from the old one, as with any simultaneous update.

Common mistakes: reading A[j][i] for A[i][j] (rows are suppliers, columns are buyers). Updating x entry by entry inside a round, so that later entries use new values. Forgetting the final demand in the update, or counting only the layers after it. Using <= where the stopping rule says <.

Practice: Orders Ripple Through the Island

The Leontief inverse: all the layers at once

The equation x = Ax + d is a linear system: (I - A)\,x = d, with I the identity matrix (1 on the diagonal, 0 elsewhere). Its solution is x = (I - A)^{-1} d = L\,d, and L = (I - A)^{-1} is the Leontief inverse, the sum of all the layers, I + A + A^2 + \cdots, in one matrix. Column j of L is the output every industry must produce, directly and down the whole chain of suppliers, to deliver one extra dollar of industry j's product to a final buyer. Its column sum, \sum_i L_{ij}, is industry j's output multiplier.

To invert a matrix by hand, extend Level 2's elimination. Gauss-Jordan elimination writes the matrix and the identity side by side, [M \mid I], and applies the same row operations to both halves until the left half is the identity; the right half is then M^{-1}. For each column c:

partial pivoting: swap in the row (from c down) with the largest entry in column c, so we never divide by 0 or by a tiny number;
divide that row by its diagonal entry, the pivot, which makes it 1;
subtract multiples of that row from every other row, above and below, to clear the rest of column c.

def inverse(M):
    n = len(M)
    W = [[float(v) for v in M[i]] + [float(i == k) for k in range(n)] for i in range(n)]   # [M | I]
    for c in range(n):
        piv = max(range(c, n), key=lambda r: abs(W[r][c]))         # partial pivoting
        if abs(W[piv][c]) < 1e-12:
            raise ValueError("the matrix cannot be inverted")
        W[c], W[piv] = W[piv], W[c]
        f = W[c][c]
        W[c] = [v / f for v in W[c]]                               # pivot becomes 1
        for r in range(n):
            if r != c:
                g = W[r][c]
                W[r] = [vr - g * vc for vr, vc in zip(W[r], W[c])]  # clear column c in row r
    return [row[n:] for row in W]

print([[round(v, 4) for v in row] for row in inverse([[0, 2], [4, 1]])])    # pivoting swaps the rows first

A = [[0.10, 0.20, 0.05], [0.15, 0.10, 0.30], [0.20, 0.25, 0.10]]
d = [50, 30, 20]
n = len(A)
I_minus_A = [[float(i == j) - A[i][j] for j in range(n)] for i in range(n)]
L = inverse(I_minus_A)
print("L =", [[round(v, 4) for v in row] for row in L])
check = [[sum(L[i][k] * I_minus_A[k][j] for k in range(n)) for j in range(n)] for i in range(n)]
print("largest error in L (I - A) = I:", max(abs(check[i][j] - (i == j)) for i in range(n) for j in range(n)) < 1e-12)
print("output multipliers", [round(sum(L[i][j] for i in range(n)), 4) for j in range(n)])
print("L d =", [round(sum(L[i][j] * d[j] for j in range(n)), 4) for i in range(n)])

[[-0.125, 0.25], [0.5, 0.0]]
L = [[1.2017, 0.3147, 0.1717], [0.3188, 1.308, 0.4537], [0.3556, 0.4333, 1.2753]]
largest error in L (I - A) = I: True
output multipliers [1.8761, 2.056, 1.9007]
L d = [72.9614, 64.2551, 56.2845]

L d is the sum of all the layers of the previous section, found without iterating. Read the tools column (index 1): a dollar of tools sold to a final buyer needs 1.31 dollars of tools in all (the dollar itself plus the tools used to make grain, transport and more tools), 0.31 of grain and 0.43 of transport, 2.06 dollars of output across the economy. Every entry of L is non-negative here, as it must be when all column sums of A are below 1: more final demand never means less output anywhere. A row sum of L answers a different question (how much industry i must make if every final demand rises by a dollar). Weighting a column by another per-dollar figure, such as jobs, emissions or water per dollar of each industry's output, gives the jobs, emissions or water behind one dollar of final demand.

Common mistakes: inverting A instead of I - A. Applying a row operation to the left half only. Clearing only the rows below the pivot (that is Gaussian elimination, which then needs back substitution; Gauss-Jordan clears above too). Summing rows when the column sums are asked for. Changing the caller's matrix instead of a copy.

Practice: Jobs Behind a Thousand Dollars

Determinants, productive economies and cost prices

Can the economy deliver at all? If the industries use up too much of each other's output, no non-negative production plan meets a positive final demand. Then I - A is singular, or its inverse has negative entries. The Hawkins-Simon condition decides it: the economy is productive exactly when every leading principal minor of I - A is positive. The k-th leading principal minor is the determinant of the top-left k \times k block.

Determinants by elimination. Reduce the matrix to an upper triangle with row operations. Subtracting a multiple of one row from another leaves the determinant unchanged; swapping two rows flips its sign; and the determinant of a triangle is the product of its diagonal. So the determinant is the product of the pivots, with one sign flip per swap, and 0 if some column has no non-zero pivot left.

When the question is "is this exactly zero?", floats are the wrong tool: the same number reached along two different routes can differ in the last bit. fractions.Fraction computes exactly. Build fractions from integers or strings, never from floats: Fraction(0.1) is the binary number nearest to 0.1, not \frac{1}{10}.

from fractions import Fraction

def det(M):
    M = [row[:] for row in M]
    n, d = len(M), 1
    for c in range(n):
        piv = next((r for r in range(c, n) if M[r][c] != 0), None)
        if piv is None:
            return 0 * d                                  # no pivot: the determinant is 0
        if piv != c:
            M[c], M[piv] = M[piv], M[c]
            d = -d                                        # a swap flips the sign
        d = d * M[c][c]
        for r in range(c + 1, n):
            g = M[r][c] / M[c][c]
            M[r] = [x - g * y for x, y in zip(M[r], M[c])]
    return d

print(det([[0, 2, 1], [3, 1, 0], [1, 1, 1]]), " (one swap)")
print(Fraction(0.1) == Fraction(1, 10), Fraction("0.1") == Fraction(1, 10))

a = [["0.7", "0.6"], ["0.45", "0.1"]]
floats = [[float(i == j) - float(a[i][j]) for j in range(2)] for i in range(2)]
exact = [[Fraction(int(i == j)) - Fraction(a[i][j]) for j in range(2)] for i in range(2)]
print("2nd minor in floats", det(floats), " exactly", det(exact))

-4.0  (one swap)
False True
2nd minor in floats 6.66133814775094e-17  exactly 0

The float minor is a tiny positive number, which would call this economy productive; exactly, it is 0, and it is not.

Physical units. When A is measured in tonnes instead of dollars, a column may sum to more than 1, and the economy can still be productive. A tonne of steel may need 1.6 tonnes of ore, as long as ore needs little steel. With ore and steel as goods 1 and 2:

from fractions import Fraction as F

def det(M):
    M = [row[:] for row in M]
    n, d = len(M), F(1)
    for c in range(n):
        piv = next((r for r in range(c, n) if M[r][c] != 0), None)
        if piv is None:
            return F(0)
        if piv != c:
            M[c], M[piv] = M[piv], M[c]
            d = -d
        d *= M[c][c]
        for r in range(c + 1, n):
            g = M[r][c] / M[c][c]
            M[r] = [x - g * y for x, y in zip(M[r], M[c])]
    return d

def minors(a):
    B = [[F(int(i == j)) - a[i][j] for j in range(len(a))] for i in range(len(a))]
    return [det([row[:k] for row in B[:k]]) for k in range(1, len(a) + 1)]

mild = [[F(1, 10), F(16, 10)], [F(5, 100), F(2, 10)]]      # column 2: 1.6 t ore and 0.2 t steel per t of steel
greedy = [[F(1, 10), F(16, 10)], [F(4, 10), F(3, 10)]]     # ore now needs 0.4 t of steel per tonne
print("mild  ", [str(m) for m in minors(mild)])
print("greedy", [str(m) for m in minors(greedy)])

mild   ['9/10', '16/25']
greedy ['9/10', '-1/100']

The first economy is productive although making steel uses 1.8 tonnes of goods per tonne; in the second, ore and steel use up so much of each other that the second minor is negative.

Cost prices: the same matrix, read the other way. In a competitive economy each price just covers the cost of a unit: the inputs it uses, at the going prices, plus the value added v_j per unit (wages, rent). Good j's inputs are its recipe, column j:

p_j = \sum_i p_i\, a_{ij} + v_j, \qquad \text{that is} \qquad p = p\,a + v

with p a row vector. Rearranged, p\,(I - a) = v; written as an ordinary system with the unknowns in a column, that is (I - a)^{\mathsf T} p = v, the transpose (rows and columns swapped). Quantities use I - a, prices use its transpose: the two are duals. For the mild economy, with value added 2 coins per tonne of ore and 5 per tonne of steel, Cramer's rule (each unknown is a determinant with its column replaced by the right-hand side, divided by the determinant of the matrix) gives:

from fractions import Fraction as F

a = [[F(1, 10), F(16, 10)], [F(5, 100), F(2, 10)]]
v = [F(2), F(5)]
B = [[F(int(i == j)) - a[i][j] for j in range(2)] for i in range(2)]
T = [[B[j][i] for j in range(2)] for i in range(2)]          # the transpose of I - a
D = T[0][0] * T[1][1] - T[0][1] * T[1][0]
p = [(v[0] * T[1][1] - T[0][1] * v[1]) / D, (T[0][0] * v[1] - v[0] * T[1][0]) / D]
print("prices", [str(x) for x in p], [round(float(x), 4) for x in p])
for j in range(2):
    print("good", j + 1, ": inputs", str(sum(p[i] * a[i][j] for i in range(2))), "+ value added", str(v[j]), "=", str(p[j]))

prices ['185/64', '385/32'] [2.8906, 12.0312]
good 1 : inputs 57/64 + value added 2 = 185/64
good 2 : inputs 225/32 + value added 5 = 385/32

A tonne of steel costs 1.6 tonnes of ore at the ore price, plus 0.2 tonnes of steel at the steel price, plus 5 coins, and the price is exactly that. In a productive economy with non-negative value added, the prices are never negative.

Common mistakes: solving (I - a)\,p = v (the rows) instead of the transpose. Deciding "exactly zero" in floats. Building a Fraction from a float. Forgetting the sign flip of a row swap. Treating a column sum above 1 in physical units as proof that the economy fails.

Practice: What Should a Crate Cost?

Many goods at once: equilibrium as a linear system

Back to exchange, now with m goods. Bisection worked for one unknown price; with many it does not, and tâtonnement can be slow. With Cobb-Douglas tastes there is an exact route, because the value of each demand is linear in the prices. Household i spends s_{ij} m_i on good j, and its income m_i = \sum_k p_k e_{ik} is linear in the prices, so the value of all demand for good j is

p_j \sum_i x_{ij} = \sum_i s_{ij} \sum_k p_k\, e_{ik} = \sum_k C_{jk}\, p_k, \qquad C_{jk} = \sum_i s_{ij}\, e_{ik}

Market j clears when demand equals the total E_j; multiplied by p_j, that is E_j\, p_j - \sum_k C_{jk}\, p_k = 0. These are m linear equations in the m prices, but they cannot pin all of them down: doubling every price still solves them, and by Walras' law any m - 1 of them imply the last. So the full matrix is singular, its determinant 0. Fix the numéraire, p_1 = 1, drop one equation, and move the p_1 terms to the right-hand side: m - 1 equations in m - 1 unknowns. Three households trade wax, flax and tar:

from fractions import Fraction as F

weights = [[2, 1, 1], [1, 3, 1], [1, 1, 2]]
owns = [[6, 0, 0], [0, 4, 2], [1, 1, 5]]
m = 3
E = [sum(e[j] for e in owns) for j in range(m)]
C = [[sum(F(w[j], sum(w)) * e[k] for w, e in zip(weights, owns)) for k in range(m)] for j in range(m)]
M = [[(E[j] if j == k else 0) - C[j][k] for k in range(m)] for j in range(m)]
print("rows of M:", [[str(v) for v in row] for row in M])
D3 = (M[0][0] * (M[1][1] * M[2][2] - M[1][2] * M[2][1]) - M[0][1] * (M[1][0] * M[2][2] - M[1][2] * M[2][0])
      + M[0][2] * (M[1][0] * M[2][1] - M[1][1] * M[2][0]))
print("determinant of the full system:", D3)

# p1 = 1; keep the equations of goods 2 and 3; unknowns p2, p3
a11, a12, r1 = M[1][1], M[1][2], -M[1][0]
a21, a22, r2 = M[2][1], M[2][2], -M[2][0]
D = a11 * a22 - a12 * a21
p = [F(1), (r1 * a22 - a12 * r2) / D, (a11 * r2 - r1 * a21) / D]
print("prices", [str(x) for x in p])
for j in range(m):
    demand = sum(F(w[j], sum(w)) * sum(pk * ek for pk, ek in zip(p, e)) / p[j] for w, e in zip(weights, owns))
    print("good", j + 1, " demand", demand, " exists", E[j])

rows of M: [['15/4', '-21/20', '-33/20'], ['-7/4', '47/20', '-49/20'], ['-2', '-13/10', '41/10']]
determinant of the full system: 0
prices ['1', '161/86', '93/86']
good 1  demand 7  exists 7
good 2  demand 5  exists 5
good 3  demand 7  exists 7

The full 3 \times 3 system has determinant exactly 0, as Walras' law predicts. After fixing the wax price and dropping the wax equation, two equations remain, and the prices they give clear all three markets, the dropped one included. With fractions the check is exact; with floats and Gaussian elimination (with partial pivoting) for larger systems, check against a tolerance instead.

Common mistakes: keeping all m equations and fixing p_1 as well, which is m equations for m - 1 unknowns. Forgetting the numéraire, which leaves a singular system. Mixing up the indices of C: row j is the good the money is spent on, column k the good whose sale earned it. Forgetting to normalise the weights.

Practice: Prices for Every Stall at Once

Production and a tax: one island in general equilibrium

The last micro model has production. One household owns the only firm on an island; coconuts are the numéraire (price 1) and the wage is w coconuts an hour. Two optimisers meet in the labour market.

The firm hires L hours, harvests Q = A\,L^{\beta} coconuts (0 < \beta < 1), and maximises its profit \pi = A L^{\beta} - w L. Like Level 3's firm it adds hours while an hour adds more than it costs, so it stops where the marginal product of labour equals the wage: by Level 4's power rule, \beta A L^{\beta - 1} = w. The profit goes to the household.

The household has T hours to split between work and leisure \ell, and Cobb-Douglas tastes u = c^{a} \ell^{1 - a}. Leisure is a good with a price: an hour of leisure is an hour not worked, and costs the wage it would have earned after tax, the net wage (1 - t)\,w. So value all T hours at the net wage and add the other income, the profit \pi and any lump-sum transfer R: that is the full income M = (1-t)\,w\,T + \pi + R. Level 2's fixed shares then give

c = a\,M, \qquad \ell = \frac{(1 - a)\,M}{(1 - t)\,w}

and, as in the Edgeworth box, the household's MRS between leisure and consumption, \frac{1 - a}{a} \cdot \frac{c}{\ell}, equals the price ratio, the net wage. Check the leisure formula against a search, for a household facing a net wage of 2 with 6 coconuts of other income:

a, T, net, other = 0.4, 16, 2.0, 6.0
best = max((k / 1000 for k in range(1, 16000)), key=lambda l: (net * (T - l) + other) ** a * l ** (1 - a))
M = net * T + other
print("search", best, " formula", round((1 - a) * M / net, 6), " consumption", round(net * (T - best) + other, 6), "=", round(a * M, 6))

search 11.4  formula 11.4  consumption 15.2 = 15.2

Equilibrium without a tax. The wage must make the hours the firm wants equal the hours the household offers. A higher wage means fewer hours wanted and, here, more offered, so bisect on the wage:

A, beta, a, T = 6, 0.6, 0.4, 16

def firm(w):                                     # hours wanted and profit at wage w
    L = (beta * A / w) ** (1 / (1 - beta))       # beta A L**(beta - 1) = w, solved for L
    return L, A * L ** beta - w * L

def offered(w, other):                           # no tax: the net wage is w
    M = w * T + other
    return T - (1 - a) * M / w

lo, hi = 1e-6, 1e6
while hi - lo > 1e-12 * hi:
    w = (lo + hi) / 2
    L, profit = firm(w)
    if L > offered(w, profit):                   # the firm wants more hours: the wage must rise
        lo = w
    else:
        hi = w
L, profit = firm(lo)
Q, M = A * L ** beta, lo * T + profit
c, leisure = a * M, T - L
print("wage", round(lo, 4), " hours", round(L, 4), " output", round(Q, 4), " profit", round(profit, 4))
print("consumption", round(c, 4), "= output", round(Q, 4), "  MRS", round((1 - a) / a * c / leisure, 4), "= wage = marginal product", round(beta * A * L ** (beta - 1), 4))

wage 1.9601  hours 4.5714  output 14.9342  profit 5.9737
consumption 14.9342 = output 14.9342   MRS 1.9601 = wage = marginal product 1.9601

Everything harvested is eaten, as it must be: the household's income is the wage bill plus the profit, which is the whole harvest. And the MRS equals the marginal product: the last hour of work produces exactly what the household thinks an hour of leisure is worth. That is efficiency, the first welfare theorem with production.

A wage tax. Now the government taxes wages at the rate t and hands the revenue back as a lump sum R. The household takes R as given (it does not see that its own extra hour raises the transfer), but in equilibrium R must equal the revenue, t\,w\,L. So in a simulation R is one more number that has to come out consistent, like the wage. The household's income barely changes, but its incentive does: it compares leisure with the net wage (1 - t)\,w, while an hour of work still produces w. That gap is the tax wedge:

\text{MRS} = (1 - t)\,w \;<\; w = \text{marginal product}

The household works too little: one more hour would produce more than the household values the hour of leisure it gives up. The cost is a deadweight loss, and the natural unit is coconuts. The consumption equivalent of a loss is the number of extra coconuts that would make the household, keeping the leisure it has, exactly as well off as before: solve (c + x)^{a}\,\ell^{1 - a} = u_0 for x. For a household that ends with 11 coconuts and 10 hours of leisure, against a utility u_0 reached with 13 coconuts and 9 hours:

a = 0.4
u0 = 13 ** a * 9 ** (1 - a)
c, leisure = 11, 10
x = (u0 / leisure ** (1 - a)) ** (1 / a) - c
print("u0", round(u0, 6), " now", round(c ** a * leisure ** (1 - a), 6), " extra coconuts", round(x, 6), " check", round((c + x) ** a * leisure ** (1 - a), 6))

u0 10.426123  now 10.388601  extra coconuts 0.099595  check 10.426123

Common mistakes: letting the household count the transfer as a reward for working. Forgetting the profit in the household's income. Valuing leisure at the gross wage when there is a tax. Measuring the consumption equivalent with the old leisure instead of the new. Forgetting that the transfer must equal the revenue in equilibrium.

Practice: A Wage Tax on Coconut Island

Time: today, tomorrow and the lifetime budget

Now the macro half, where choices are spread over periods t = 0, 1, 2, \ldots. Period 0 is today: its money is not discounted. Level 2's present value turns a whole stream into today's money, and a household that can save and borrow at the rate r carries a balance a_t from period to period. A common timing: in period t it receives its income y_t, spends c_t, and whatever is left (negative if it borrowed) earns interest until the next period,

a_{t+1} = (1 + r)\,(a_t + y_t - c_t)

Follow the balance and compare it with present values:

r, a0 = 0.04, 15.0
incomes = [0, 30, 50, 50, 10]
spend = [20, 25, 25, 25, 25]
a, balances = a0, [a0]
for y, c in zip(incomes, spend):
    a = (1 + r) * (a + y - c)
    balances.append(a)
print("balances", [round(b, 4) for b in balances])
T = len(incomes)
pv_y = sum(y / (1 + r) ** t for t, y in enumerate(incomes))
pv_c = sum(c / (1 + r) ** t for t, c in enumerate(spend))
print("PV of income", round(pv_y, 4), " PV of spending", round(pv_c, 4))
print("final balance", round(balances[-1], 4), "= (1 + r)**T * (a0 + PV income - PV spending) =", round((1 + r) ** T * (a0 + pv_y - pv_c), 4))

balances [15.0, -5.2, -0.208, 25.7837, 52.815, 39.3276]
PV of income 128.0718  PV of spending 110.7474
final balance 39.3276 = (1 + r)**T * (a0 + PV income - PV spending) = 39.3276

The household is in debt after periods 0 and 1, then saves. The last line is the lifetime budget: whatever the path in between, the final balance is the starting wealth plus the present value of income minus the present value of spending, carried forward to period T. Ending with exactly nothing, a_T = 0, means spending has the same present value as lifetime wealth W = a_0 + \sum_t y_t/(1 + r)^t. A constant spending level c has present value c times the sum of the discount factors, \sum_{t=0}^{T-1} (1 + r)^{-t}.

Two different numbers discount the future, and they must not be mixed up. The market turns a coin next period into \frac{1}{1 + r} of a coin today: that is a price, set by the interest rate. A person's patience is a discount factor \beta: one unit of happiness next period is worth \beta units now. Budgets use 1 + r; lifetime utility uses \beta.

Sums like 1 + q + q^2 + \cdots + q^{T-1} have the closed form \frac{1 - q^T}{1 - q} (Level 2's geometric series), which divides by zero when q = 1, that is when r = 0 or \beta = 1. A loop has no such case:

def sum_powers(q, T):                    # 1 + q + q**2 + ... + q**(T - 1)
    total, term = 0.0, 1.0
    for t in range(T):
        total += term
        term *= q
    return total

for q in (1 / 1.04, 0.9, 1):
    formula = round((1 - q ** 10) / (1 - q), 6) if q != 1 else "division by zero"
    print(round(q, 6), " loop", round(sum_powers(q, 10), 6), " formula", formula)

0.961538  loop 8.435332  formula 8.435332
0.9  loop 6.513216  formula 6.513216
1  loop 10.0  formula division by zero

Common mistakes: discounting period 0. Adding interest before the period's income and spending when the stated timing says the leftover earns interest. Dividing wealth by T to get a level spending, which is right only when r = 0. Discounting a budget with \beta, or utility with 1 + r. Using the closed form when q may be 1.

Practice: Spreading a Boatbuilder's Earnings

The Euler equation: moving one coin cannot help

How should spending be spread over time? Take log utility: a period with consumption c is worth \ln c, and a plan is worth \sum_t \beta^t \ln c_t. The logarithm rises ever more slowly, so the first coins of a period matter most, which pushes towards smooth plans. Two periods, wealth W today, and whatever is not spent today earns interest:

\max_{c_0}\; \ln c_0 + \beta \ln c_1, \qquad c_1 = (1 + r)(W - c_0)

A search over c_0:

import math
W, r, beta = 100, 0.05, 0.9

def U(c0):
    return math.log(c0) + beta * math.log((1 + r) * (W - c0))

c0 = max((k / 1000 for k in range(1, 100000)), key=U)
c1 = (1 + r) * (W - c0)
print("best c0", c0, " c1", round(c1, 4), " c1 / c0", round(c1 / c0, 4), " beta (1 + r)", round(beta * (1 + r), 6))
for plan in (c0, 60.0):
    eps = 0.01                                   # move one cent from period 0 to period 1
    print("from c0 =", plan, " moving a cent changes U by", round(U(plan - eps) - U(plan), 8))

best c0 52.632  c1 49.7364  c1 / c0 0.945  beta (1 + r) 0.945
from c0 = 52.632  moving a cent changes U by -3e-08
from c0 = 60.0  moving a cent changes U by 5.829e-05

The best plan grows by the factor \beta(1 + r) from one period to the next. Here is why. Move one small coin \varepsilon from period 0 to period 1. Today you lose about \varepsilon \cdot \frac{1}{c_0} of utility (the slope of \ln at c_0 is 1/c_0). Tomorrow you have (1 + r)\varepsilon more, worth \beta (1 + r) \varepsilon \cdot \frac{1}{c_1} in today's terms. At the best plan the move can neither help nor hurt, so the two are equal:

\frac{1}{c_t} = \beta\,(1 + r)\,\frac{1}{c_{t+1}} \quad\Longleftrightarrow\quad c_{t+1} = \beta\,(1 + r)\,c_t

This is the Euler equation for log utility. At the best plan, moving a cent changes U only in the eighth decimal (by a second-order amount); from the worse plan c_0 = 60 it gains at once. The Euler equation fixes the shape of a plan with any number of periods: each period's consumption is the previous one times g = \beta(1 + r). Impatience beats interest when g < 1, and consumption falls; interest wins when g > 1. The lifetime budget of the previous section then fixes the level: one equation in the one unknown c_0.

Saving in a growing economy. The same logic decides how much a whole economy should save. Take Level 4's Solow economy, y = A k^{\alpha}, with full depreciation (all capital wears out within a period, so next period's capital is this period's saving, k' = s\,y) and log utility of consumption c = (1 - s)\,y. Simulate many periods for each constant saving rate s and score each with \sum_t \beta^t \ln c_t:

import math
A, alpha, beta = 2, 0.35, 0.8

def lifetime(s, k0, T=300):
    k, total = k0, 0.0
    for t in range(T):
        y = A * k ** alpha
        total += beta ** t * math.log((1 - s) * y)
        k = s * y
    return total

for k0 in (0.1, 5.0):
    best = max((s / 1000 for s in range(1, 1000)), key=lambda s: lifetime(s, k0))
    print("start", k0, " best constant saving rate", best, "  alpha * beta =", round(alpha * beta, 6))

start 0.1  best constant saving rate 0.28   alpha * beta = 0.28
start 5.0  best constant saving rate 0.28   alpha * beta = 0.28

Whatever the starting capital, the best rate is \alpha\beta. A fuller argument (guessing that the value of capital is a constant plus a multiple of \ln k, and checking the guess in the Bellman equation of a later section) shows that no changing plan does better: with log utility and full depreciation, the optimal saving rate is exactly \alpha\beta, at every capital level. Solow's assumption of a fixed saving rate happens to be optimal here. With capital that lasts longer, it is not.

Common mistakes: writing the Euler equation upside down, c_{t+1} = c_t / (\beta(1 + r)). Using \beta in the budget or 1 + r in the utility. Concluding that a flat plan is best: it is only when \beta(1 + r) = 1. Forgetting that the lifetime budget is still needed to find the level.

Practice: Patience Against Interest

Backward induction: value functions with a last day

Many dynamic choices are made step by step, each depending on where the last one left you. A value function V_t(x) answers: if at time t you are in state x (the jars left, the cash in hand), how well can you do from now on, choosing well? On the day after the last, nothing is left to gain, so V_T(x) = 0. On any earlier day, the best you can do is the best choice now plus the discounted value of where it leaves you, the Bellman equation:

V_t(x) = \max_{q}\; \Big[\, \text{reward}_t(q) + \beta\, V_{t+1}(\text{next state}) \,\Big]

So fill the tables backwards, from day T - 1 to day 0, storing for every state both the value and the best choice (the policy). This is the dynamic programming of the Algorithms path: every subproblem (a day and a state) is solved once, from smaller ones that are already known.

A beekeeper has 6 jars of honey to sell over three market days. Prices differ by day, and the more she sells on one day, the lower the price she gets: selling q jars on day t earns q\,(\text{price}_t - q). Unsold jars are worth nothing after the last day, and she discounts with \beta = 0.9:

N, prices, beta = 6, [10, 14, 9], 0.9
T = len(prices)

def revenue(t, q):
    return q * (prices[t] - q)

V = [0.0] * (N + 1)                         # V_T: jars left after the last day are worth 0
policy = [None] * T
for t in range(T - 1, -1, -1):
    newV, sell = [0.0] * (N + 1), [0] * (N + 1)
    for x in range(N + 1):
        vals = [revenue(t, q) + beta * V[x - q] for q in range(x + 1)]
        top = max(vals)
        q = next(q for q, v in enumerate(vals) if v >= top - 1e-9)    # the smallest best choice
        newV[x], sell[x] = vals[q], q
    V, policy[t] = newV, sell
    print("day", t, " V", [round(v, 2) for v in V], " sell", sell)

x, plan = N, []
for t in range(T):                          # follow the stored choices forwards
    plan.append(policy[t][x])
    x -= policy[t][x]
print("plan", plan, " worth", round(V[N], 4))

day 2  V [0.0, 8.0, 14.0, 18.0, 20.0, 20.0, 20.0]  sell [0, 1, 2, 3, 4, 4, 4]
day 1  V [0.0, 13.0, 24.0, 33.0, 40.2, 47.2, 52.6]  sell [0, 1, 2, 3, 3, 4, 4]
day 0  V [0.0, 11.7, 21.6, 30.6, 38.7, 45.7, 52.18]  sell [0, 0, 0, 1, 1, 2, 2]
plan [2, 3, 1]  worth 52.18

The last day's table is simple: sell until one more jar would lower the price by more than it brings. With 5 or 6 jars left on day 2, selling 4 and selling 5 both earn exactly 20, and the tie rule below picks 4. Each earlier day weighs selling now against the discounted value of keeping jars. The plan is read forwards: start with 6 jars on day 0, take that day's stored choice, move to the state it leaves, and look up the next day's own table.

Ties. Two choices can be equally good in real arithmetic and differ in the last bits of their floats, because their sums were added in a different order. So never break ties with ==. Find the best value first, then take the first choice (here the smallest) whose value is within a stated tolerance of it, as the next(...) line does: "the smallest choice within 1e-9 of the best". When a problem asks for a whole plan with a tie rule, applying the rule day by day while following the stored choices forwards gives the earliest-smallest plan.

Common mistakes: filling the tables forwards. Building the plan from the day-0 table only. Forgetting the choice of 0. Breaking ties with == or with >, which keeps the largest choice. Updating V in place while the same day's other states still need the old values.

Practice: A Wheel of Cheese for the Winter

For ever: value function iteration

Without a last day there is no table to start from. A household that plans for ever, or a fishery that should last for ever, faces the same problem every period, so its value function does not depend on t:

V(x) = \max_{q}\; \Big[\, \text{reward}(q) + \beta\, V(\text{next state}) \,\Big]

V appears on both sides. Value function iteration (VFI) treats the right-hand side as a machine that turns a guess into a better guess: start from V_0 = 0 everywhere, compute V_1 from V_0 for every state at once, then V_2 from V_1, and so on. V_k is exactly the value of a problem with k periods left, so VFI is backward induction run until the tables stop changing.

Why does it stop changing? Take two guesses that differ by at most D in every state. After one update their difference is at most \beta D: everything they disagree about sits behind a factor \beta. The update is a contraction, the largest change shrinks at least by the factor \beta each time, and the guesses converge to the one solution from any start. Stop on the largest absolute change over all states, \max_x |V_{k+1}(x) - V_k(x)|, and count the update that met the tolerance. The policy is read off that last update.

A fishing village manages a stock of x fish, in hundreds, from 0 to 12. Each season it catches c (from 0 to x), worth \sqrt{c}, and the fish left behind breed: next season's stock is min(12, 3 * (x - c) // 2). With \beta = 0.9:

import math
N, beta, tol = 12, 0.9, 1e-6

def regrow(left):
    return min(N, 3 * left // 2)

V = [0.0] * (N + 1)
updates, changes = 0, []
while True:
    newV, policy = [0.0] * (N + 1), [0] * (N + 1)          # a whole new table from the old one
    for x in range(N + 1):
        vals = [math.sqrt(c) + beta * V[regrow(x - c)] for c in range(x + 1)]
        top = max(vals)
        c = next(c for c, v in enumerate(vals) if v >= top - 1e-9)
        newV[x], policy[x] = vals[c], c
    updates += 1
    change = max(abs(u - v) for u, v in zip(newV, V))
    changes.append(change)
    V = newV
    if change < tol:
        break
print("updates", updates, " catch for each stock", policy)
print("first changes", [round(c, 4) for c in changes[:5]], " later ratios", [round(changes[k + 1] / changes[k], 4) for k in (20, 40, 60)])
x, path = 2, [2]
for season in range(8):
    x = regrow(x - policy[x])
    path.append(x)
print("stock from 2:", path)

updates 138  catch for each stock [0, 1, 0, 0, 0, 1, 1, 1, 2, 1, 2, 3, 4]
first changes [3.4641, 1.6854, 1.4851, 1.3365, 1.2029]  later ratios [0.9, 0.9, 0.9]
stock from 2: [2, 3, 4, 6, 7, 9, 12, 12, 12]

The changes shrink by exactly \beta = 0.9 per update once the policy has settled, and the loop stops after the update whose change first falls below the tolerance. From a small stock the village fishes lightly or not at all and lets the stock grow; at the top it takes the surplus the stock would not keep. The policy comes from the last update, since earlier updates belong to shorter horizons.

Two details matter in code. First, never update the table in place: entries computed in this update would be read as if they were old ones, and the count of updates (which problems check) would change. Second, a reward table that does not change between updates (here \sqrt{c}; elsewhere \ln of consumption) is cheaper computed once; skip the choices that are not allowed (consumption 0 under log utility) instead of taking log(0).

Common mistakes: updating V in place. Stopping on the change in one state, or on a signed change. Counting one update too few (the update that met the tolerance counts). Returning the policy of the update before the last. Mixing up grid numbers that start at 1 with list positions that start at 0.

Practice: Until the Values Stop Moving, The Island Chooses Its Own Saving Rate

Random futures: expected values in the Bellman equation

So far the future was known. When next period's income or offer is random, with values y_i and probabilities p_i, the Bellman equation uses the expected value of where a choice leads: the probability-weighted average over what might happen,

V_t(m) = \max_{s}\; \Big[\, u(m - s) + \beta \sum_i p_i\, V_{t+1}(\text{next state if } y_i) \,\Big]

The choice is made before the luck is known, and each later choice may depend on how the luck turned out, which is why a value function is computed for every state, not one fixed plan.

Risk and saving. With log utility, risk on its own makes a household save more. A farmer with 10 coins of cash this season, a gross return R = 1.1 on saving, \beta = 0.95, and next season's harvest worth either 6 for sure, or 1 or 11 with equal chances (also 6 on average):

import math
m, R, beta = 10, 1.1, 0.95
for name, ys, ps in (("safe", [6], [1]), ("risky", [1, 11], [0.5, 0.5])):
    def value(s):
        return math.log(m - s) + beta * sum(p * math.log(R * s + y) for y, p in zip(ys, ps))
    best = max((k / 100 for k in range(0, 1000)), key=value)
    print(name, " best saving", best, " value", round(value(best), 4))

s = 2
print("E[ln(Rs + y)]", round(0.5 * math.log(R * s + 1) + 0.5 * math.log(R * s + 11), 4), " < ln(Rs + E[y])", round(math.log(R * s + 6), 4))
print("E[1/(Rs + y)]", round(0.5 / (R * s + 1) + 0.5 / (R * s + 11), 4), " > 1/(Rs + E[y])", round(1 / (R * s + 6), 4))

safe  best saving 2.07  value 4.0785
risky  best saving 3.29  value 3.9045
E[ln(Rs + y)] 1.8717  < ln(Rs + E[y]) 2.1041
E[1/(Rs + y)] 0.1941  > 1/(Rs + E[y]) 0.122

The average harvest is the same, yet the risky farmer saves more than twice as much. The last two lines show why. The logarithm bends, so the average of \ln is below \ln of the average: risk costs utility. More importantly, the marginal value of a saved coin, the average of 1/(\text{next cash}), is higher with risk, because a bad harvest makes a coin very valuable: this is precautionary saving. Taking the average income first and then the value would miss it entirely. A borrowing limit (saving at least 0) adds a corner: a household that expects a large income would like to borrow against it, cannot, and saves 0.

Waiting for a better offer. Some problems have one number at their heart. A boat owner gets one offer a week: 300, 400, 500 or 600 coins, with probabilities 0.1, 0.4, 0.3 and 0.2, independently each week. Selling ends the story with the offer in hand. Waiting costs 20 coins of upkeep this week and brings a fresh offer next week, discounted by \beta = 0.95. Let h be the value of waiting. Each offer is worth the better of selling and waiting, \max(o, h), and h is the upkeep plus the discounted average of next week's offer values:

h = -20 + \beta \sum_i p_i \max(o_i, h)

h is on both sides. Solve it by repeated substitution: guess h, compute the right-hand side, use it as the next guess. Each round shrinks the error by at least \beta (the same contraction as VFI), so the guesses converge:

offers, probs, cost, beta = [300, 400, 500, 600], [0.1, 0.4, 0.3, 0.2], 20, 0.95
h, rounds, changes = 0.0, 0, []
while True:
    new = -cost + beta * sum(p * max(o, h) for o, p in zip(offers, probs))
    rounds += 1
    changes.append(abs(new - h))
    h = new
    if changes[-1] < 1e-10:
        break
print("rounds", rounds, " value of waiting", round(h, 6), " sells at", [o for o in offers if o >= h])
print("first changes", [round(c, 2) for c in changes[:4]], " ratio late on", round(changes[-1] / changes[-2], 4))
print("check: -20 + 0.95 * (0.5 h + 0.3 * 500 + 0.2 * 600) =", round(-20 + 0.95 * (0.5 * h + 150 + 120), 6))

rounds 37  value of waiting 450.47619  sells at [500, 600]
first changes [417.0, 17.57, 8.35, 3.97]  ratio late on 0.4752
check: -20 + 0.95 * (0.5 h + 0.3 * 500 + 0.2 * 600) = 450.47619

The owner sells for 500 or 600 and turns down 300 and 400. The cutoff h plays the part of a reservation value: accept every offer at or above it. When an offer is exactly worth h the owner is indifferent, and a problem must say what happens then (usually: accept). Once you know which offers are accepted, the equation is linear and can be checked by hand, as the last line does. In job search the same idea gives a reservation wage; there a job accepted pays its wage every period for ever, so its value is a geometric series, w + \beta w + \beta^2 w + \cdots = \frac{w}{1 - \beta}.

Common mistakes: averaging the income before applying the utility. Letting the household borrow when the problem forbids it, or save in the last period. Stopping substitution after a fixed number of rounds (patient searchers need hundreds). Valuing a job at one period's wage instead of the whole stream. Dividing by zero when no offer is ever accepted.

Practice: A Rainy-Day Jar for the Ferry Pilot, Holding Out for a Better Offer

The cobweb: planting on last year's price

Not every model assumes foresight. Farmers must plant months before they sell, and in the cobweb model they plant as if this year's price will be last year's. The harvest of season t is set by P_{t-1}; once picked it must be sold, so the price P_t is whatever makes buyers take exactly the harvest:

a - b\,P_t = c + d\,P_{t-1}

Pig farmers, with demand Q = 90 - 3P and supply Q = 6 + 2.4\,P_{t-1} (tonnes, coins per kilo):

a, b, c, d = 90, 3, 6, 2.4
p_star = (a - c) / (b + d)                     # the price that repeats itself
P, prices = 12.0, [12.0]
for season in range(8):
    harvest = c + d * P                         # planted on last season's price
    P = (a - harvest) / b                       # the price at which buyers take the harvest
    prices.append(P)
gaps = [p - p_star for p in prices]
print("p* =", round(p_star, 4))
print("prices", [round(p, 3) for p in prices])
print("gaps  ", [round(g, 3) for g in gaps])
print("ratios", [round(gaps[t + 1] / gaps[t], 4) for t in range(4)], " -d/b =", round(-d / b, 6))

p* = 15.5556
prices [12.0, 18.4, 13.28, 17.376, 14.099, 16.721, 14.623, 16.301, 14.959]
gaps   [-3.556, 2.844, -2.276, 1.82, -1.456, 1.165, -0.932, 0.746, -0.597]
ratios [-0.8, -0.8, -0.8, -0.8]  -d/b = -0.8

A cheap year makes farmers plant little; the small harvest fetches a high price; they plant a lot, and the price falls: a zigzag. The equilibrium price p^* is the one that repeats itself, P_t = P_{t-1} = p^*, so a - b\,p^* = c + d\,p^*. Subtract that equation from the season's equation and the constants cancel:

b\,(P_t - p^*) = -d\,(P_{t-1} - p^*) \quad\Longrightarrow\quad P_t - p^* = -\frac{d}{b}\,(P_{t-1} - p^*)

The gap is multiplied by -d/b every season, a first-order difference equation like Level 2's output catching up, but with a negative factor, hence the zigzag. If d < b the swings shrink, if d = b they repeat for ever at the same size, if d > b they grow. That verdict belongs to the market, not to the path: a market started exactly at p^* never moves, even when any other start would swing further and further away. Press Run to see the cobweb that gives the model its name:

a, b, c, d = 90, 3, 6, 2.4
P, xs, ys = 12.0, [], []
for season in range(10):
    q = c + d * P
    xs += [q, q]
    ys += [P, (a - q) / b]
    P = (a - q) / b
qs = [k for k in range(20, 70)]
plot(qs, [(a - q) / b for q in qs], label="demand", xlabel="tonnes", ylabel="price")
plot(qs, [(q - c) / d for q in qs], label="supply (last season's price)")
plot(xs, ys, label="the path")

Aside (optional): naive or rational expectations? The cobweb's farmers never notice that their forecast is wrong in the same direction every other year. John Muth's answer, which grew into rational expectations, was that people use the model itself to forecast, so forecasts are right on average: here they would expect p^*, and the price would sit there. Supporters see this as a discipline (a model should not assume people make the same mistake for ever) and it underlies most modern macroeconomics. Critics reply that it asks a great deal of knowledge and computing from ordinary people, and point to experiments and surveys in which forecasts do lag behind, adapt slowly, or follow rules of thumb. Models with learning or several kinds of forecasters sit in between. Problems here state the expectation rule they use.

Common mistakes: using this season's price on the supply side, which is an ordinary market that clears at once. Judging the verdict from a path that starts at p^*. Comparing d with b the wrong way round: steep supply (large d) over-reacts. Returning T prices when the path from season 0 to season T has T + 1.

Practice: Cabbages Planted on Last Year's Price

Crowds: thresholds and cascades

Some outcomes come from people watching each other. Two simple models show how a crowd can tip all at once, or settle on a choice most of its members privately doubt.

Who reacts to what, and when. A crowd simulation in rounds must say what each person sees: usually the state at the end of the previous round, so everyone in a round reacts to the same numbers. Then a round first collects all its changes, and only afterwards applies them, as in VFI. Updating while scanning lets a change be seen in its own round, and the result then depends on the order of the list. A row of six lamps: lamp 0 is on, and in each round an off lamp switches on if the lamp to its left was on:

on = [True] + [False] * 5
rounds = 0
while True:
    new = [i for i in range(1, 6) if not on[i] and on[i - 1]]    # decided from the old state
    if not new:
        break
    for i in new:
        on[i] = True
    rounds += 1
print("simultaneous:", rounds, "rounds")

on = [True] + [False] * 5
for i in range(1, 6):                                           # in place: each change is seen at once
    if not on[i] and on[i - 1]:
        on[i] = True
print("in place: everything in 1 round", on)

simultaneous: 5 rounds
in place: everything in 1 round [True, True, True, True, True, True]

In a threshold model (Mark Granovetter's) person i acts, joins a protest or withdraws savings, once at least k_i others have acted. A crowd whose thresholds run 0, 1, 2, 3, ... tips one person per round until everybody has acted; take out the single person with threshold 1, and the crowd stays almost still. Two crowds with nearly the same thresholds can end in opposite places.

Information cascades. Customers choose one of two stalls, A and B, one at a time. Each has a private signal, a hint pointing to A or B; all signals are equally reliable and independent. Everyone sees the earlier choices, not the earlier signals. With equally reliable signals, beliefs depend only on a count: the A signals a customer knows of minus the B signals. She adds her own signal (+1 for A, -1 for B) to the public lead and goes with the sign; at 0, a stated rule decides (say, she follows her own signal).

What can later customers learn from her choice? Only what it depends on. Try both possible signals: if they lead to different choices, her choice reveals her signal and the public lead moves by 1; if both lead to the same choice, it reveals nothing:

def choose(lead, signal):                     # signal +1 (A) or -1 (B); a tie follows the own signal
    total = lead + signal
    if total == 0:
        return "A" if signal > 0 else "B"
    return "A" if total > 0 else "B"

for lead in range(-3, 4):
    with_a, with_b = choose(lead, 1), choose(lead, -1)
    print("lead", lead, " signal A ->", with_a, " signal B ->", with_b, " reveals the signal:", with_a != with_b)

lead -3  signal A -> B  signal B -> B  reveals the signal: False
lead -2  signal A -> B  signal B -> B  reveals the signal: False
lead -1  signal A -> A  signal B -> B  reveals the signal: True
lead 0  signal A -> A  signal B -> B  reveals the signal: True
lead 1  signal A -> A  signal B -> B  reveals the signal: True
lead 2  signal A -> A  signal B -> A  reveals the signal: False
lead 3  signal A -> A  signal B -> A  reveals the signal: False

With a lead of -1, 0 or 1 every choice reveals a signal. At \pm 2 both signals give the same choice: the customer herds, her choice carries no information, the lead freezes, and the next customer faces exactly the same situation. A cascade has begun, and it lasts for ever unless someone's signal becomes public some other way. For the signals B, B, A, A, A, A:

customer	signal	public lead before	choice	revealed?	lead after
1	B	0	B	yes	-1
2	B	-1	B	yes	-2
3	A	-2	B	no (herds)	-2
4 to 6	A	-2	B	no (herd)	-2

Four of six signals pointed to A, and everybody went to B. Two early signals decided the whole queue, because later customers rationally ignored their own information, and so never passed it on.

Common mistakes: counting every earlier choice as a signal (choices made in a cascade reveal nothing). Breaking a tie towards the crowd when the rule says the own signal. Applying round changes while still scanning the round. Counting the final round in which nothing happens.

Practice: Queue at the Savings Bank, Following the Queue at the Noodle Stalls

Agent-based models: a market of individuals

An agent-based model simulates individual traders following simple rules, and measures what the market as a whole achieves. There is no equation to solve; the model is the procedure, so the work is to follow it to the letter. Level 5's rules for seeded simulation become strict:

create one generator, rng = random.Random(seed), and draw every random number from it;
draw in the stated order (all values, then all costs; shuffle buyers, then sellers);
shuffle lists of agent numbers and look their values up, so that removing an agent is easy;
when the procedure says the remaining agents go back into increasing order, rebuild the list with sorted(...).

rng.shuffle uses the same random numbers for any list of the same length, so shuffling numbers or shuffling values gives the same order of agents; but shuffling a list in a different order (sorted against unsorted) gives a different market:

import random
names = ["ana", "ben", "cy", "dee", "eli"]
order = list(range(5))
random.Random(3).shuffle(names)
random.Random(3).shuffle(order)
print(names, order)

['ana', 'cy', 'dee', 'eli', 'ben'] [0, 2, 3, 4, 1]

A labour market with random meetings: five workers with skills and five vacancies with requirements, drawn from one generator. Each round, the workers still looking and the vacancies still open are shuffled and paired in order; a pair matches when the skill is at least the requirement, and both leave:

import random
rng = random.Random(35)
skills = [rng.randint(1, 10) for _ in range(5)]          # workers first,
needs = [rng.randint(1, 10) for _ in range(5)]           # then vacancies
print("skills", skills, " needs", needs)
workers, jobs, matches = list(range(5)), list(range(5)), []
for rnd in range(3):
    rng.shuffle(workers)                                 # workers first, then vacancies
    rng.shuffle(jobs)
    hired, filled = set(), set()
    for w, j in zip(workers, jobs):                      # zip pairs in order and stops at the shorter list
        if skills[w] >= needs[j]:
            matches.append((w, j))
            hired.add(w)
            filled.add(j)
    workers = sorted(w for w in workers if w not in hired)
    jobs = sorted(j for j in jobs if j not in filled)
    print("round", rnd + 1, " matches so far", matches, " still looking", workers, " still open", jobs)

# the most matches possible: fill the easiest vacancies with the least skilled workers who qualify
ws, best, i = sorted(skills), 0, 0
for need in sorted(needs):
    while i < len(ws) and ws[i] < need:
        i += 1
    if i < len(ws):
        best += 1
        i += 1
print("matches", len(matches), " most possible", best, " efficiency", round(len(matches) / best, 4))

skills [9, 6, 3, 6, 3]  needs [5, 7, 5, 10, 1]
round 1  matches so far [(0, 2), (4, 4)]  still looking [1, 2, 3]  still open [0, 1, 3]
round 2  matches so far [(0, 2), (4, 4), (3, 0)]  still looking [1, 2]  still open [1, 3]
round 3  matches so far [(0, 2), (4, 4), (3, 0)]  still looking [1, 2]  still open [1, 3]
matches 3  most possible 4  efficiency 0.75

Random meetings waste chances. In round 1 the most skilled worker (skill 9) met a vacancy needing only 5, which a worker with skill 6 could have filled; after that nobody left can fill the vacancy needing 7, and from round 3 on the market is stuck at 3 matches. Pairing on sorted lists would have made 4. The model's efficiency compares what the market achieved with the best possible outcome, computed separately (here by a greedy match on sorted lists; for a market of buyers and sellers, by pairing the highest values with the lowest costs). Markets with richer rules, posted prices, haggling or memory, are compared the same way: realised over maximal.

Common mistakes: creating a generator inside the loop, or using the module's random functions. Drawing in a different order than stated. Leaving the remaining agents in shuffled order when the procedure says increasing order. Pairing values directly, so that removing agents shifts the lists. Dividing by a best outcome of 0.

Practice: Random Pairs in the Grain Pit

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

Checklist

You can now:

Value an endowment at market prices, compute Cobb-Douglas demands and excess demands, and use Walras' law as a check.
Fix a numéraire and find a clearing price by bisection or by algebra, and run a tâtonnement loop in its stated order.
Compute marginal rates of substitution in an Edgeworth box, say who gains by buying good 1, trace the contract curve, and compare allocations with a tolerance under the first welfare theorem.
Read an input-output matrix by columns, sum the layers d + Ad + A^2 d + \cdots by iteration, and invert I - A by Gauss-Jordan elimination with partial pivoting to get output multipliers.
Compute determinants by elimination, exactly with Fraction when a zero matters, test the Hawkins-Simon condition, and solve the dual price equation with the transpose.
Turn Cobb-Douglas market clearing for many goods into a linear system, drop one equation and fix the numéraire.
Solve a production economy with a firm, a household and a tax, explain the tax wedge, and measure a loss as a consumption equivalent.
Follow a balance with stated timing, use the lifetime budget, keep \beta and 1 + r apart, and derive the Euler equation by moving one coin.
Solve finite-horizon problems by backward induction and infinite-horizon ones by value function iteration, with expected values over a random future, a tolerance-based tie rule, and fixed points found by repeated substitution.
Simulate the cobweb and judge its stability from -d/b, track what a crowd can infer in an information cascade, and run seeded agent-based markets measured against their best possible outcome.