Functions, modules, and program structure

Alfred Galichon (NYU)

'math+econ+code' masterclass series: fundamentals of research in python

With python code examples

© 2018–2026 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.

If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/

Learning objectives

• Define functions with positional, keyword, default, and variadic (*args, **kwargs) parameters, and know when each is the right tool.

• Recognize Python's scope rules, the LEGB chain, and the role of global and nonlocal.

• Use lambdas, map, filter, and functools.partial sparingly and idiomatically.

• Organize code across files: write a module, import it cleanly, and understand how packages, environments, and if __name__ == "__main__" fit together.

• Build an operator, a function that returns a function, and locate its fixed point. This is the structure of the Solow growth model here, and of the Bellman operator in the dp series.

• Report and check a solver's convergence status rather than silently reading a number out of a routine that may not have converged.

• Recognize when a class or a @dataclass earns its keep.

References

[M] McKinney, W. (2022). Python for Data Analysis (3rd ed.), Chapter 3. https://wesmckinney.com/book/.

[Py] The Python Tutorial, §§4, 6, 9. https://docs.python.org/3/tutorial/.

[QE] Sargent, T. J. and Stachurski, J. QuantEcon Python Lectures. https://python-programming.quantecon.org/.

[S] Solow, R. M. (1956). "A Contribution to the Theory of Economic Growth." Quarterly Journal of Economics 70(1), 65–94: the model of §9.

[BS] Barro, R. J. and Sala-i-Martin, X. (2004). Economic Growth (2nd ed.). MIT Press, Chapter 1: the convergence rate computed in §9 and Exercise 4.

1. Motivation

fd01 gave us expressions and fd02 gave us containers. Neither gave us a way to name a computation and reuse it, which is what separates a script from a piece of research software.

This lecture is about that, and it has a specific target. The 'math+econ+code' method is to write a generic solver first, and then applications that reuse it: not to re-derive the algorithm inside each model. That requires functions to be first-class objects: passed as arguments (a solver takes the function whose root it seeks), returned from other functions (an operator is built from parameters and returns a map), and packaged into modules that other notebooks import.

Section 8 writes our first such module, mec_numerical, with a Newton solver and a bisection solver. Section 9 then uses it on the Solow growth model, which is this lecture's worked example and its economics. The pattern established there, parameters in, operator out, fixed point of the operator, convergence rate checked against theory, is the same one that returns in the dp series with the Bellman operator, where the fixed point is the value function.

Remark on the duality lens. This lecture contains no optimization problem, hence no multiplier and no dual variable, and it would be dishonest to manufacture one. What it contributes to the series' spine is the representation: an operator as a first-class object. The dual objects of the later lectures, value functions, potentials, prices, are all computed as fixed points of operators built exactly this way.

2. Defining functions

A function is a named, reusable computation. Three things make it work: the def keyword, the parameters named in the signature, and what is returned.

def discount_factor(rate, periods):
    """Return the price today of one unit delivered after `periods` periods."""
    return (1 + rate) ** (-periods)

discount_factor(0.05, 10)

0.6139132535407591

The triple-quoted string at the top is a docstring. It serves three purposes: it documents the function, it appears under help(discount_factor) and discount_factor? in Jupyter, and tools such as Sphinx and pdoc extract it to build API documentation. Every function you intend anyone else to read should have one, and "anyone else" includes you, six months from now.

Note what this particular function is: the price $p_t = (1+r)^{-t}$ of fd01 §11, the first dual object of the series, now given a name and a signature. Naming it is the whole point of this lecture.

help(discount_factor)

Help on function discount_factor in module __main__:

discount_factor(rate, periods)
    Return the price today of one unit delivered after `periods` periods.

A function with no return returns None: the same sentinel met in fd01. Functions that exist for their side effects (printing, mutating a container, writing a file) return it implicitly.

3. Default and keyword arguments

Parameters can carry default values. Callers may then pass them positionally, by keyword, or omit them.

def newton(f, df, x0, tol=1e-10, maxiter=100):
    """Solve f(x) = 0 by Newton's method. Returns (root, iterations, converged)."""
    x = x0
    for k in range(maxiter):
        fx = f(x)
        if abs(fx) < tol:
            return x, k, True
        x = x - fx / df(x)
    return x, maxiter, False        # exhausted the budget without converging

# all positional
newton(lambda x: x ** 2 - 2, lambda x: 2 * x, 1.0)

(1.4142135623746899, 4, True)

Return the convergence status, and check it. The third return value is not decoration. A solver that has run out of iterations still returns a number, and that number is not a root; a caller who unpacks only the first element will carry it forward as though it were. This is one of the standard ways a numerical result silently goes wrong, and it is why the house rule is to check the status before reading the solution: a rule that applies equally to scipy.optimize in fd09 and to linprog in the lp series, both of which report status codes that are easy to ignore.

The cell below gives the solver a budget it cannot meet, and shows the difference between the two habits.

f, df = lambda x: x ** 2 - 2, lambda x: 2 * x

root, k, converged = newton(f, df, 1.0, maxiter=2)      # deliberately too few
print(f"returned root = {root:.10f}   after {k} iterations   converged = {converged}")
print(f"residual |f(root)| = {abs(f(root)):.2e}   <- not a root")

if not converged:
    print("\nstatus checked: this result must not be used.")

returned root = 1.4166666667   after 2 iterations   converged = False
residual |f(root)| = 6.94e-03   <- not a root

status checked: this result must not be used.

Keyword-only parameters. Anything after a bare * in the signature must be passed by name. Use this when an argument's meaning is not obvious from its position: verbose=True is the classic case, and so is any boolean flag.

def bisect(f, a, b, *, tol=1e-10, maxiter=200, verbose=False):
    """Bisection on a sign-change interval. Returns (root, iterations, converged)."""
    fa = f(a)
    if fa * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")
    for k in range(maxiter):
        m = 0.5 * (a + b)
        fm = f(m)
        if verbose:
            print(f"  k={k:3d}  m={m:.8f}  f(m)={fm:+.2e}")
        if abs(fm) < tol or 0.5 * (b - a) < tol:
            return m, k, True
        if fa * fm < 0:
            b = m
        else:
            a, fa = m, fm
    return m, maxiter, False

bisect(f, 1.0, 2.0, verbose=False)

(1.4142135623842478, 28, True)

Note the guard: bisect raises rather than proceeding when the bracket does not contain a sign change. Failing loudly at the point of the mistake is worth far more than returning a plausible number computed from invalid assumptions: the same argument as for KeyError over a silent default in fd02 §5.

Defaults are evaluated once, at definition time. If a default is a mutable object, every call that omits the argument shares it: the trap of fd02 §7. In a solver signature this matters concretely: a default like history=[] would accumulate the iterates of every call ever made. Use None as the sentinel and build the container inside.

4. *args and **kwargs

A starred parameter collects any number of positional arguments into a tuple; a double-starred parameter collects keyword arguments into a dict. The names args and kwargs are convention: the * and ** are what matter. This is the collecting side of the operators met in fd02 §11.

def trace(label, *args, **kwargs):
    print(f"[{label}] positional: {args}")
    print(f"[{label}] keyword:    {kwargs}")

trace("demo", 1, 2, 3, name="newton", tol=1e-8)

[demo] positional: (1, 2, 3)
[demo] keyword:    {'name': 'newton', 'tol': 1e-08}

On the call side, * and ** unpack a sequence or a dict into arguments. This is how numerical configuration is forwarded through a chain of functions without restating it at every level: one dict of settings, defined once, passed down.

config = {"tol": 1e-6, "maxiter": 20}

root, k, converged = newton(f, df, 1.0, **config)
print(f"root = {root:.10f}  iterations = {k}  converged = {converged}")

root = 1.4142135624  iterations = 4  converged = True

5. Scope: the LEGB rule

When Python meets a name inside a function it searches four scopes, in order:

1. Local: the function's own scope;
2. Enclosing: the scope of any function wrapping it;
3. Global: the module's top level;
4. Built-in: print, len, range, and the rest.

The first match wins.

r_global = 0.05          # global

def pv_annuity(c, T):
    return c * (1 - (1 + r_global) ** -T) / r_global      # reads the global

pv_annuity(100, 5)

432.9476670630823

Reading versus writing. A function reads names from enclosing scopes without ceremony. To write to a name in an outer scope you must declare it: global for module level, nonlocal for an enclosing function. Without the declaration, an assignment creates a new local of the same name, and the outer one is untouched.

counter = 0

def bump():
    global counter          # without this, the next line would create a local
    counter += 1

bump(); bump(); bump()
counter

3

Treat global as a code smell. A function that mutates module-level state is hard to test, because its result depends on history rather than on its arguments; and it is hard to reason about, because any caller anywhere can change what it does. The cleaner pattern is to return the new value and let the caller rebind the name.

The example above also explains something about pv_annuity: it silently depends on r_global. Change that global anywhere in the notebook and the function's answer changes, with nothing in its signature to warn you. A function whose behaviour is fully determined by its arguments is testable; one that reads globals is not. This is exactly the property fd04 will rely on when we write a test suite.

6. Lambdas, map, and filter

A lambda is an anonymous, single-expression function. Use one when a small function is needed as an argument and naming it would only clutter the page: as with the f and df passed to newton above.

gdp_2022 = [("USA", 25.46), ("China", 17.96), ("Japan", 4.23), ("Germany", 4.07)]
sorted(gdp_2022, key=lambda pair: pair[1], reverse=True)

[('USA', 25.46), ('China', 17.96), ('Japan', 4.23), ('Germany', 4.07)]

map(f, xs) and filter(pred, xs) apply a function or a predicate across an iterable. Both return iterators, not lists, the laziness of fd02 §9, so wrap them in list(...) to materialize. In modern Python a comprehension usually reads better:

xs = range(1, 11)

via_map = list(map(lambda x: x ** 2, filter(lambda x: x % 2 == 0, xs)))
via_comprehension = [x ** 2 for x in xs if x % 2 == 0]

print(via_map)
print(via_comprehension)
print("identical:", via_map == via_comprehension)

[4, 16, 36, 64, 100]
[4, 16, 36, 64, 100]
identical: True

When the callable already has a name: str.upper, math.sqrt: map(name, xs) is perfectly good style. Reserve lambdas for ad-hoc keys and one-line arguments; anything longer than a single short expression deserves a def, a name, and a docstring.

7. Closures and functools.partial

A function defined inside another function closes over the enclosing scope: it remembers the names of the function that created it, even after that function has returned. This is a closure, and this construction is central to the lecture.

def make_discounter(rate):
    """Return the price function t -> (1 + rate)^(-t) for a fixed rate."""
    def price(t):
        return (1 + rate) ** (-t)
    return price

p_lo = make_discounter(0.02)        # two different price functions,
p_hi = make_discounter(0.10)        # from one factory

print(f"price of a dollar in 10 years at  2%: {p_lo(10):.6f}")
print(f"price of a dollar in 10 years at 10%: {p_hi(10):.6f}")

price of a dollar in 10 years at  2%: 0.820348
price of a dollar in 10 years at 10%: 0.385543

Read what just happened economically. make_discounter takes a rate and returns the entire price system of fd01 §11 as a single object. The parameters are separated from the map they induce, and the map can then be passed around, applied, or handed to a solver. That separation, parameters in, operator out, is the shape of every structural model in this masterclass series:

$$ \text{primitives} \;\longmapsto\; \text{an operator} \;\longmapsto\; \text{its fixed point} \;=\; \text{the equilibrium.} $$

In §9 the primitives are $(s, \alpha, \delta)$ and the operator is the Solow accumulation map. In dp01 the primitives are $(\Phi, P, \beta)$ and the operator is the Bellman operator, whose fixed point is the value function.

functools.partial does a related job: it fixes some arguments of an existing function rather than building a new one from scratch.

from functools import partial

def pv(c, r, T):
    return c * (1 - (1 + r) ** -T) / r

pv_at_5pct = partial(pv, r=0.05)

print(f"partial : {pv_at_5pct(c=100, T=5):.6f}")
print(f"direct  : {pv(c=100, r=0.05, T=5):.6f}")
print("identical:", pv_at_5pct(c=100, T=5) == pv(c=100, r=0.05, T=5))

partial : 432.947667
direct  : 432.947667
identical: True

Closures can also carry state, through nonlocal. The counter below wraps a function and records how many times it is called: the standard way to measure the cost of a solver in function evaluations, which is the honest currency when each evaluation is expensive. You will build one in Exercise 3.

def counted(func):
    """Wrap func, counting calls. Returns (wrapped, get_count)."""
    calls = 0
    def wrapped(x):
        nonlocal calls          # rebind the enclosing 'calls', do not shadow it
        calls += 1
        return func(x)
    return wrapped, lambda: calls

f_counted, count = counted(f)
newton(f_counted, df, 1.0)
print(f"Newton used {count()} evaluations of f")

Newton used 5 evaluations of f

8. Modules, packages, and environments

A module is a .py file. Importing it makes its names available:

import math                  # then math.exp, math.sqrt, ...
from math import exp, sqrt   # exp and sqrt directly
import numpy as np           # the standard alias

Three conventions you will meet constantly:

• prefer import x to from x import *; the latter pollutes the namespace and hides where a name came from;
• library aliases are fixed by custom: numpy as np, pandas as pd, matplotlib.pyplot as plt, tensorflow as tf, and readers expect them;
• imports go at the top of the file in three blank-line-separated groups: standard library, third-party, local.

Writing your own module. A module is just a file, and the %%writefile cell magic writes the rest of the cell to disk. We use it to assemble the small numerical toolbox that §9 and the exercises will import.

import os
os.makedirs("generated", exist_ok=True)

%%writefile generated/mec_numerical.py
"""mec_numerical: a small root-finding toolbox for the fd series.

Every solver returns (root, iterations, converged). Callers are expected to
check `converged` before using `root`.
"""


def newton(f, df, x0, tol=1e-10, maxiter=100):
    """Solve f(x) = 0 by Newton's method from x0."""
    x = x0
    for k in range(maxiter):
        fx = f(x)
        if abs(fx) < tol:
            return x, k, True
        x = x - fx / df(x)
    return x, maxiter, False


def bisect(f, a, b, tol=1e-10, maxiter=200):
    """Solve f(x) = 0 by bisection on a bracket [a, b] with a sign change."""
    fa = f(a)
    if fa * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")
    for k in range(maxiter):
        m = 0.5 * (a + b)
        fm = f(m)
        if abs(fm) < tol or 0.5 * (b - a) < tol:
            return m, k, True
        if fa * fm < 0:
            b = m
        else:
            a, fa = m, fm
    return m, maxiter, False


def fixed_point(operator, x0, tol=1e-12, maxiter=10_000):
    """Iterate x <- operator(x) from x0 until successive iterates agree."""
    x = x0
    for k in range(maxiter):
        x_next = operator(x)
        if abs(x_next - x) < tol:
            return x_next, k, True
        x = x_next
    return x, maxiter, False


if __name__ == "__main__":
    # smoke test: runs only when this file is executed as a script
    root, k, ok = newton(lambda x: x * x - 2, lambda x: 2 * x, 1.0)
    print(f"sqrt(2) = {root} in {k} iterations (converged={ok})")

Writing generated/mec_numerical.py

The file now sits in generated/. Adding that directory to the import path lets us import it like any other module.

import sys
if "generated" not in sys.path:
    sys.path.insert(0, "generated")

import mec_numerical

root, k, converged = mec_numerical.newton(f, df, 1.0)
print(f"sqrt(2) = {root:.12f}  ({k} iterations, converged={converged})")
print(f"module docstring: {mec_numerical.__doc__.splitlines()[0]}")

sqrt(2) = 1.414213562375  (4 iterations, converged=True)
module docstring: mec_numerical: a small root-finding toolbox for the fd series.

The if __name__ == "__main__": idiom. When Python imports a file it executes it top to bottom. Any code at module level therefore runs on every import, which is rarely what you want for tests or demonstrations. The guard runs that code only when the file is executed as a script: python mec_numerical.py, and not when it is imported. Note that the smoke test above did not print when we imported the module, which is the guard doing its job.

Two related points:

• Packages. A package is a directory of modules containing an __init__.py (possibly empty), which lets you write from mypkg.numerical import newton. We will not need one until the code outgrows a single file.
• importlib.reload. Re-importing a module already imported in the same kernel does not re-read the file; Python caches it. Use importlib.reload(mec_numerical) to pick up an edit, or, cleaner, and the only version that is reproducible, restart the kernel and run all.

Environments. Every project gets its own. The two standard recipes:

# conda, recommended for scientific work
conda create -n mec python=3.12 numpy scipy pandas matplotlib jupyterlab
conda activate mec
conda env export --from-history > environment.yml

# plain pip + venv
python -m venv .venv
source .venv/bin/activate            # Windows: .venv\Scripts\activate
pip install numpy scipy pandas matplotlib jupyterlab
pip freeze > requirements.txt

The exported environment.yml or requirements.txt is what makes a result reproducible on someone else's machine, and it is the first item on the reproducibility checklist we assemble in fd07. Conda resolves the compiled dependencies underneath NumPy and SciPy more cleanly than pip alone; either is fine, but be in an environment: never install into the system Python.

9. Worked example: the Solow model as an operator

We now have everything needed for the lecture's economics. In the Solow (1956) model, output per worker is $f(k) = k^\alpha$ with $0 < \alpha < 1$, a constant fraction $s$ of it is saved, and capital depreciates at rate $\delta$. Capital per worker then evolves according to

$$ k_{t+1} \;=\; T(k_t) \;:=\; s\,k_t^{\alpha} \;+\; (1-\delta)\,k_t . \tag{9.1} $$

$T$ is an operator: given the primitives $(s, \alpha, \delta)$ it is a map from capital today to capital tomorrow. In code it is exactly the closure of §7: parameters in, function out.

A steady state is a fixed point of $T$, a level $k^\ast$ with $T(k^\ast) = k^\ast$; equivalently, investment just covers depreciation:

$$ s\,k^{\ast\alpha} = \delta k^\ast \qquad\Longleftrightarrow\qquad k^\ast = \left(\frac{s}{\delta}\right)^{\frac{1}{1-\alpha}} . \tag{9.2} $$

We will find $k^\ast$ three independent ways: by iterating the operator, by solving $g(k) = s k^\alpha - \delta k = 0$ with Newton, and from the closed form (9.2), and then check something sharper than the level: the rate at which the economy approaches it.

def make_solow_step(s, alpha, delta):
    """Return the Solow accumulation operator T(k) for given primitives."""
    def step(k):
        return s * k ** alpha + (1 - delta) * k
    return step

s, alpha, delta = 0.2, 0.3, 0.05
T = make_solow_step(s, alpha, delta)

k_star_closed = (s / delta) ** (1 / (1 - alpha))
print(f"primitives: s = {s}, alpha = {alpha}, delta = {delta}")
print(f"closed form  k* = (s/delta)^(1/(1-alpha)) = {k_star_closed:.12f}")
print(f"is it a fixed point?  T(k*) - k* = {T(k_star_closed) - k_star_closed:+.2e}")

primitives: s = 0.2, alpha = 0.3, delta = 0.05
closed form  k* = (s/delta)^(1/(1-alpha)) = 7.245789314111
is it a fixed point?  T(k*) - k* = +0.00e+00

Route 1: iterate the operator. Start the economy poor, at $k_0 = 0.1$, and apply $T$ until successive iterates stop moving. This is fixed_point from the module we just wrote.

Route 2: Newton on $g(k) = sk^{\alpha} - \delta k$, whose positive root is $k^\ast$ by (9.2). One warning about the starting point, and it is a live instance of the caveat in fd01 §12: Newton's guarantees are local, and $g$ is nearly flat near $k = 1$, so a step from there overshoots badly. We start at $k_0 = 5$ instead, and check below what the other choice would have done.

k_iter, n_iter, converged_iter = mec_numerical.fixed_point(T, x0=0.1, tol=1e-12)

print(f"iteration : k* = {k_iter:.12f}   after {n_iter} periods   converged = {converged_iter}")
assert converged_iter, "the iteration did not converge"

# route 2: k* is the root of g(k) = s k^alpha - delta k, solved by Newton
def g(k):
    return s * k ** alpha - delta * k

def dg(k):
    return s * alpha * k ** (alpha - 1) - delta

k_newton, n_newton, converged_newton = mec_numerical.newton(g, dg, x0=5.0)
print(f"Newton    : k* = {k_newton:.12f}   after {n_newton} iterations  converged = {converged_newton}")
assert converged_newton, "Newton did not converge"

# route 3: the closed form
print(f"closed form: k* = {k_star_closed:.12f}")

gap = max(abs(k_iter - k_star_closed), abs(k_newton - k_star_closed))
tol = 1e-8
print(f"\nmax gap across the three routes = {gap:.2e}   (tolerance {tol:.0e})")
assert gap < tol, "the three routes disagree"
print("check passed: three routes, one steady state.")

iteration : k* = 7.245789314084   after 747 periods   converged = True
Newton    : k* = 7.245789314111   after 4 iterations  converged = True
closed form: k* = 7.245789314111

max gap across the three routes = 2.69e-11   (tolerance 1e-08)
check passed: three routes, one steady state.

What a bad start would have done. The cell below takes a single Newton step from $k_0 = 1$.

k0 = 1.0
lands_at = k0 - g(k0) / dg(k0)
print(f"g(1) = {g(k0):+.6f},   g'(1) = {dg(k0):+.6f}   <- nearly flat")
print(f"one Newton step from k0 = 1 lands at k = {lands_at:.4f}")
print(f"\nbut capital cannot be negative, and k^alpha is not real there:")
print(f"  (-14.0) ** 0.3 = {(-14.0) ** 0.3}")

g(1) = +0.150000,   g'(1) = +0.010000   <- nearly flat
one Newton step from k0 = 1 lands at k = -14.0000

but capital cannot be negative, and k^alpha is not real there:
  (-14.0) ** 0.3 = (1.297349820271451+1.7856488371481518j)

At $k=1$ the marginal product of capital almost exactly offsets depreciation, so $g'(1) \approx 0.01$; dividing by it throws the iterate to $k \approx -14$. That is outside the model's domain, capital cannot be negative, and Python does not object: it evaluates $(-14)^{0.3}$ as a complex number rather than raising. The iteration then wanders through the complex plane and, in this instance, returns to the correct real answer carrying a zero imaginary part. That is worse than failing outright, because the output looks very nearly normal.

Three lessons, all of which recur. Newton is fast but not safe. Bisection on a bracket is slow but cannot leave the interval, which is why Exercise 2 pairs the two and why fd09 reaches for Brent's method, which combines them. And a solver knows nothing about the domain of your model: keeping $k > 0$ is economics, and it is your job to impose it.

Three routes agree, but look at the second column: Newton needed 5 iterations and the economic iteration needed several hundred periods. That is not a defect of the code. It is the model's own statement about how fast an economy converges, and it deserves to be checked against theory rather than merely observed.

The convergence rate. Near the steady state, linearizing (9.1) gives $k_{t+1} - k^\ast \approx T'(k^\ast)(k_t - k^\ast)$, so the error contracts geometrically at rate $T'(k^\ast)$. Differentiating and using $s k^{\ast\alpha-1} = \delta$ from (9.2),

$$ T'(k^\ast) \;=\; s\alpha k^{\ast\alpha-1} + 1 - \delta \;=\; \alpha\delta + 1 - \delta \;=\; 1 - \delta(1-\alpha). \tag{9.3} $$

A sharp, parameter-free prediction: the error ratio should converge to $1 - \delta(1-\alpha)$, here $1 - 0.05 \times 0.7 = 0.965$. Note that this is linear convergence, a constant factor per step, in contrast to the quadratic convergence verified for Newton in fd01 §12, where the number of correct digits doubled. As in that check, the measurement is destroyed by cancellation once the error approaches machine precision, so we stop the table while the ratio is still measurable.

predicted = 1 - delta * (1 - alpha)

# re-run the iteration, keeping the whole path
path = [0.1]
for _ in range(1200):
    path.append(T(path[-1]))

print(f"predicted rate T'(k*) = 1 - delta(1-alpha) = {predicted:.9f}\n")
print(f"{'t':>5}{'k_t':>15}{'growth of k':>14}{'error':>11}{'e_(t+1)/e_t':>14}{'gap':>10}")
print("-" * 69)
for t in (0, 1, 5, 20, 100, 200, 300, 400, 500):
    e_t = path[t] - k_star_closed
    if abs(e_t) < 1e-8:               # below this the ratio is rounding noise
        break
    ratio = (path[t + 1] - k_star_closed) / e_t
    growth = (path[t + 1] - path[t]) / path[t]
    print(f"{t:>5}{path[t]:>15.8f}{100 * growth:>12.3g}%{e_t:>11.2e}"
          f"{ratio:>14.9f}{abs(ratio - predicted):>10.1e}")

# the asymptotic ratio must match the theory
e_400 = path[400] - k_star_closed
ratio_400 = (path[401] - k_star_closed) / e_400
gap, tol = abs(ratio_400 - predicted), 1e-6
print("-" * 69)
print(f"at t = 400: error = {e_400:.1e}, ratio = {ratio_400:.9f}, "
      f"predicted = {predicted:.9f}, gap = {gap:.1e}   (tolerance {tol:.0e})")
assert gap < tol, "the observed convergence rate does not match the theory"
print("check passed: the economy converges at the rate the theory predicts.")

predicted rate T'(k*) = 1 - delta(1-alpha) = 0.965000000

    t            k_t   growth of k      error   e_(t+1)/e_t       gap
---------------------------------------------------------------------
    0     0.10000000        95.2%  -7.15e+00   0.986672229   2.2e-02
    1     0.19523745        57.8%  -7.05e+00   0.984007514   1.9e-02
    5     0.70735480        20.5%  -6.54e+00   0.977838531   1.3e-02
   20     2.88856489        4.52%  -4.36e+00   0.970047522   5.0e-03
  100     6.96194930       0.142%  -2.84e-01   0.965210348   2.1e-04
  200     7.23769067     0.00392%  -8.10e-03   0.965005872   5.9e-06
  300     7.24555958    0.000111%  -2.30e-04   0.965000166   1.7e-07
  400     7.24578280    3.15e-06%  -6.52e-06   0.965000005   4.7e-09
  500     7.24578913    8.93e-08%  -1.85e-07   0.964999999   7.9e-10
---------------------------------------------------------------------
at t = 400: error = -6.5e-06, ratio = 0.965000005, predicted = 0.965000000, gap = 4.7e-09   (tolerance 1e-06)
check passed: the economy converges at the rate the theory predicts.

Reading the result

The early rows and the late rows say different things, and both are economic.

The growth column is the prediction that made this model famous. At $t = 0$ the economy is very poor, capital is scarce, its marginal product is enormous, and capital grows by 95% in a single period. By $t = 100$ growth is a fraction of a percent, and it falls to zero as the steady state is approached. Poor economies grow fast and rich ones grow slowly: not by assumption, but because $f'(k) = \alpha k^{\alpha-1}$ declines as $k$ rises. Diminishing returns are the convergence mechanism.

The ratio column measures something different, and the two should not be confused. It reports how fast the remaining gap to the steady state closes, and it approaches $0.965$ from above: gap-closing is slowest at the start, around 1.3% per period, and accelerates to the asymptotic 3.5%. This is not in tension with the paragraph above; it is a reminder that (9.3) is a local statement. At $k_0 = 0.1$ the measured ratio is a secant slope across an enormous interval, not the derivative at the fixed point: indeed $T'(0.1) \approx 1.25 > 1$, so the map is locally expanding out there, and the fixed point attracts only in the large. An asymptotic rate has to be checked asymptotically, which is why the assertion below reads the ratio at $t = 400$ rather than at $t = 0$, and why Exercise 4 must iterate for thousands of periods in the slowly-converging cases.

The asymptotic rate itself is the headline. With $\delta = 0.05$ and $\alpha = 0.3$ the gap to the steady state closes by only $3.5\%$ per period, a half-life of

$$ \frac{\ln 2}{-\ln(1 - \delta(1-\alpha))} \;=\; \frac{\ln 2}{0.0356} \;\approx\; 19 \text{ periods.} $$

Read as years, an economy halfway from its steady state takes about two decades to close half the remaining distance, and the eight-digit agreement we demanded took some seven hundred years. Convergence in the neoclassical growth model is slow, and the empirical literature that grew out of this observation, the "2% per year" convergence rate of Barro and Sala-i-Martin, is essentially the measurement of (9.3) in cross-country data. Our $3.5\%$ is the right order of magnitude, and the discrepancy is informative: matching the observed $2\%$ requires a capital share $\alpha$ substantially above the $0.3$ of national accounts, which is the standard argument for a broad notion of capital that includes human capital.

Note finally what made this check possible. We did not merely observe that the iteration converged; we predicted the rate from (9.3) and confirmed it to six decimals. A number can be right by accident, but a rate matching a formula derived independently is evidence that the operator implements the model. That is the verification habit of this series, applied to an economic model rather than to arithmetic.

10. A first look at classes

A class bundles state with behaviour. Every class has an initializer, __init__; calling the class creates an instance and runs the initializer on it. Methods take the instance as their first argument, conventionally named self.

The right moment to use a class is when you find yourself passing the same tuple of arguments through a family of related functions. Wrapping them names the bundle and lets you write obj.method(...) in place of f(obj.a, obj.b, obj.c, ...).

class Bond:
    """Plain-vanilla coupon bond: face value, coupon per period, maturity."""

    def __init__(self, face, coupon, maturity):
        self.face = face
        self.coupon = coupon
        self.maturity = maturity

    def cashflows(self):
        """The stream [c, c, ..., c + F], as in fd01 section 10."""
        flows = [self.coupon] * self.maturity
        flows[-1] += self.face
        return flows

    def pv(self, r):
        """Present value at a constant rate r."""
        return sum(c / (1 + r) ** (t + 1) for t, c in enumerate(self.cashflows()))

    def ytm(self, price, r0=0.05):
        """Yield to maturity: the rate at which pv(r) equals `price`."""
        flows = self.cashflows()
        def excess(r):
            return sum(c / (1 + r) ** (t + 1) for t, c in enumerate(flows)) - price
        def d_excess(r):
            return -sum((t + 1) * c / (1 + r) ** (t + 2) for t, c in enumerate(flows))
        root, _, converged = mec_numerical.newton(excess, d_excess, r0)
        if not converged:
            raise RuntimeError("the yield-to-maturity solve did not converge")
        return root

    def __repr__(self):
        return (f"Bond(face={self.face}, coupon={self.coupon}, "
                f"maturity={self.maturity})")


bond = Bond(face=1000, coupon=100, maturity=5)
print(bond)
print(f"cashflows      : {bond.cashflows()}")
print(f"PV at 5%       : {bond.pv(0.05):.6f}")

Bond(face=1000, coupon=100, maturity=5)
cashflows      : [100, 100, 100, 100, 1100]
PV at 5%       : 1216.473834

Note that cashflows() builds [self.coupon] * self.maturity: a list of immutable floats, so the repetition operator is safe here. Had the elements been mutable, this is exactly the aliasing bug of fd02, Exercise 5.

Note also how ytm is written: it defines the excess-value function and its derivative as closures over flows, hands them to the solver we wrote in §8, and, following §3, raises rather than returning a number the solver disowns.

print(f"PV at 5%              : {bond.pv(0.05):.4f}")
for price in (1000.0, 950.0, 1216.473834):
    y = bond.ytm(price)
    print(f"priced at {price:>10.2f} -> YTM = {y:>8.4%}   (check: PV at that rate = {bond.pv(y):.4f})")

PV at 5%              : 1216.4738
priced at    1000.00 -> YTM = 10.0000%   (check: PV at that rate = 1000.0000)
priced at     950.00 -> YTM = 11.3653%   (check: PV at that rate = 950.0000)
priced at    1216.47 -> YTM =  5.0000%   (check: PV at that rate = 1216.4738)

The third line is the verification: pricing the bond at its own yield must return the price we started from, and it does. It also reconciles this lecture with fd01 §10, where the same bond at $r = 5\%$ was worth $1216.4738$: quoted at that price, its yield to maturity is exactly $5\%$, as it must be.

The economics in the first two lines is the standard inverse relation: at par the yield equals the coupon rate of $10\%$, and a lower price implies a higher yield. Price and yield are two encodings of the same object, and ytm is the change of coordinates between them.

11. @dataclass for the record case

When a class is mostly a bundle of data, an __init__ that stores its arguments and a __repr__ that lists them, the dataclasses decorator writes the boilerplate.

There is a natural use for one right here. Our solvers return the bare tuple (root, iterations, converged), and every caller has to remember that order; unpacking it wrongly is silent. A small record type fixes that.

from dataclasses import dataclass

@dataclass
class SolverResult:
    """The outcome of a root-finding call."""
    root: float
    iterations: int
    converged: bool

    def value(self):
        """Return the root, refusing if the solver did not converge."""
        if not self.converged:
            raise RuntimeError(f"solver did not converge in {self.iterations} iterations")
        return self.root


good = SolverResult(*mec_numerical.newton(f, df, 1.0))
bad = SolverResult(*mec_numerical.newton(f, df, 1.0, maxiter=2))

print(good)
print(f"good.value() = {good.value():.12f}\n")
print(bad)
try:
    bad.value()
except RuntimeError as err:
    print("bad.value() -> RuntimeError:", err)

SolverResult(root=1.4142135623746899, iterations=4, converged=True)
good.value() = 1.414213562375

SolverResult(root=1.4166666666666667, iterations=2, converged=False)
bad.value() -> RuntimeError: solver did not converge in 2 iterations

The annotations root: float are type hints. Python does not enforce them at runtime, but mypy and modern editors use them to catch errors before the code runs: we install mypy in fd04. For a dataclass they are not optional: they are how the decorator learns what the fields are.

The design point is worth more than the syntax. SolverResult.value() makes it impossible to read a root out of a failed solve by accident: the unchecked path raises instead of returning a plausible number. Turning a convention that callers are asked to remember into a structure that enforces itself is the most reliable form of verification there is, because it does not depend on anyone remembering to check.

12. Summary

• A function is a named computation, and in Python functions are values: they can be passed to other functions (a solver receives the function whose root it seeks), returned from them (§7), and stored in modules that other notebooks import (§8). That is what makes "write the generic solver once, then reuse it" possible, and it is the organizing rule of this series.

• Check the convergence status before reading the answer. A solver out of iterations still returns a number. Our solvers return (root, iterations, converged), and §11 turns that convention into a SolverResult whose value() refuses to hand over a root the solver disowns. The same discipline applies to scipy.optimize in fd09 and to linprog in the lp series.

• A closure separates the primitives of a model from the operator they induce: parameters in, function out. make_solow_step(s, alpha, delta) returns the map $T$, and the equilibrium is its fixed point. The Bellman operator of dp01 has exactly this shape, and its fixed point is the value function.

• The economics is §9. Three independent routes, iteration, Newton, closed form, agree on the steady state $k^\ast = (s/\delta)^{1/(1-\alpha)}$, and the rate of convergence matches the predicted $1 - \delta(1-\alpha)$ to six decimals. That rate is the result: capital converges at $3.5\%$ per period, a half-life of about 19 years, and closing eight digits took seven centuries. Diminishing returns make convergence slow, and the gap between this $3.5\%$ and the $2\%$ measured across countries is the standard argument for a broad, human-capital-inclusive notion of $k$.

13. Exercises

Write your answer in the cell below each prompt. Worked solutions are in §15, at the end of this notebook.

Exercise 1: IRR as a function. In fd01 you computed an internal rate of return inline. Now package it: write

def irr(c_t, r0=0.05, tol=1e-12, maxiter=100):
    ...

returning the rate at which $\sum_t c_t (1+r)^{-t} = 0$, using mec_numerical.newton. Follow §3: return the convergence status alongside the rate, or raise if the solve fails.

Test on c_t = [-1000, 200, 300, 400, 500], and verify by an independent route: mec_numerical.bisect on a bracket: reporting the gap and your tolerance. Then check the behaviour of your function on a stream with no root, such as [-1000, 100, 100, 100]: what does it do, and what should it do?

# your answer here

Exercise 2: Configuration forwarding with `kwargs.** Writecompare_solvers(f, df, a, b, kwargs)that solves $f = 0$ both withmec_numerical.newton(started at the midpoint of $[a,b]$) and withmec_numerical.bisecton $[a,b]$, forwardingkwargs:tol,maxiter`: to both, and returns a dict reporting each root, each iteration count, and the gap between them.

Apply it to $f(x) = x^2 - 2$ on $[1, 2]$ at tol=1e-12, and to the Solow $g(k) = s k^\alpha - \delta k$ on $[1, 20]$. Which solver takes more iterations, and why does the comparison reverse nothing about which you should prefer for a problem where $f$ is expensive to evaluate?

# your answer here

Exercise 3: A closure that counts, and the true cost of a solver. Using nonlocal as in §7, write counted(func) returning a wrapped function and a way to read the call count. Use it to measure how many evaluations of $f$ and, for Newton, of $f'$ each solver needs to reach tol=1e-12 on $x^2 - 2$ over $[1,2]$.

Iteration counts are the wrong currency when an evaluation is expensive: each Newton step costs one $f$ and one $f'$, while each bisection step costs one $f$. Report evaluations, not iterations, and state which method is cheaper on this problem. Then say how your answer would change if $f'$ were unavailable and had to be approximated by a finite difference.

# your answer here

Exercise 4: The Solow convergence rate (proof, then check). Prove (9.3): that the Solow operator (9.1) satisfies

$$ T'(k^\ast) = 1 - \delta(1-\alpha) $$

at the steady state, stating where you use the steady-state condition (9.2). Deduce the half-life $\ln 2 / [-\ln T'(k^\ast)]$, and show that the rate is independent of the saving rate $s$: explain, in a sentence, why $s$ moves the level of the steady state but not the speed of convergence to it.

Then verify numerically: for each $(\alpha, \delta)$ in $\{0.3, 0.5, 0.7\} \times \{0.02, 0.05, 0.10\}$, iterate the operator, measure the asymptotic error ratio, and compare it with the prediction, reporting the worst gap across the nine cases and your tolerance. Which combination gives a half-life closest to the 2% per year found in the cross-country data?

# your answer here

Exercise 5: A @dataclass for a portfolio. Define a dataclass Position with fields ticker: str, quantity: int, price: float, and a method value() returning quantity * price. Then define Portfolio holding a list of positions, with total_value() and a weights() method returning each ticker's share of the total.

Build a three-position example and report the total two ways, by summing value() over positions, and by a comprehension over the raw fields, checking that they agree. Verify also that the weights sum to one, to a stated tolerance.

Finally, one design question: Portfolio holds a list of positions rather than inheriting from list. Given §7 of fd02, what would go wrong if two Portfolio objects were built from the same list?

# your answer here

14. Further directions

You can now write functions, close over parameters to build operators, package code into modules, and bundle state with behaviour in classes.

fd04 picks the module of §8 straight back up and makes it something a colleague could depend on: it profiles a Solow loop, compiles it with Numba, and builds a pytest suite around the convergence rate derived in §9 and the $s$-independence of Exercise 4. Its Exercise 3 returns to the Newton solver of §3 and finds a bug in it that this lecture's checks did not catch.

fd05 then steps away from the language to the workflow around it: Git, GitHub, branching, and the research compendium that makes a result something another economist can clone and re-run. The module written in §8, and the tests written around it in fd04, are the first artifacts of this series that belong under version control.

Data comes next, fd06 for tables and fd07 for acquiring them, and then the scientific stack proper. In fd08 the Solow path of §9 becomes a NumPy array computed without a Python loop, the mec_numerical solvers give way to scipy.optimize in fd09, and the operator-and-fixed-point structure of §9 reappears throughout the dp series, where $T$ is the Bellman operator, its fixed point is the value function, and the value function is the multiplier on a mass-balance constraint.

Save your work, restart the kernel, and run all cells top-to-bottom before you move on.

15. Solutions to the exercises

Reference solutions. As in the lecture, every numerical answer is checked against a second route with a stated tolerance.

Solution to Exercise 1: IRR as a function

def irr(c_t, r0=0.05, tol=1e-12, maxiter=100):
    """Internal rate of return of the stream c_t, dated t = 0, 1, ...."""
    def g(r):
        return sum(c / (1 + r) ** t for t, c in enumerate(c_t))

    def dg(r):
        return sum(-t * c / (1 + r) ** (t + 1) for t, c in enumerate(c_t))

    root, iterations, converged = mec_numerical.newton(g, dg, r0, tol=tol, maxiter=maxiter)
    if not converged:
        raise RuntimeError(f"IRR solve did not converge in {iterations} iterations")
    return root, iterations


c_t = [-1000, 200, 300, 400, 500]
rate, iterations = irr(c_t)
print(f"IRR = {rate:.12f}  ({100 * rate:.4f}% per period, {iterations} iterations)")

IRR = 0.128257269002  (12.8257% per period, 5 iterations)

# independent route: bisection, which uses no derivative
def g(r):
    return sum(c / (1 + r) ** t for t, c in enumerate(c_t))

rate_bisect, k_bisect, ok = mec_numerical.bisect(g, 0.0, 1.0, tol=1e-14)
print(f"bisection IRR = {rate_bisect:.12f}  ({k_bisect} iterations, converged={ok})")

gap, tol = abs(rate - rate_bisect), 1e-9
print(f"|Newton - bisection| = {gap:.2e}   (tolerance {tol:.0e})")
assert gap < tol, "the two solvers disagree"
print(f"residual |g(IRR)| = {abs(g(rate)):.2e}")
print("check passed.")

bisection IRR = 0.128257269002  (46 iterations, converged=True)
|Newton - bisection| = 1.69e-15   (tolerance 1e-09)
residual |g(IRR)| = 8.53e-14
check passed.

# a stream with no root: every discounted sum stays negative
hopeless = [-1000, 100, 100, 100]
print("value of the hopeless stream at several rates:")
for r in (0.0, 0.05, 0.5, 5.0):
    print(f"  r = {r:>5.2f} -> {sum(c / (1 + r) ** t for t, c in enumerate(hopeless)):>10.2f}")

try:
    irr(hopeless)
except (RuntimeError, OverflowError, ZeroDivisionError) as err:
    print(f"\nirr(hopeless) raised {type(err).__name__}: {err}")

value of the hopeless stream at several rates:
  r =  0.00 ->    -700.00
  r =  0.05 ->    -727.68
  r =  0.50 ->    -859.26
  r =  5.00 ->    -980.09

irr(hopeless) raised OverflowError: (34, 'Result too large')

What it does, and what it should do. The stream repays $300$ in total against an outlay of $1000$, so its value is negative at every rate: no IRR exists. Newton, hunting for a crossing that is not there, drives $r$ upward without bound, and the run above ends not in a tidy non-convergence report but in an OverflowError raised by $(1+r)^t$, long before the iteration budget is exhausted.

That is worth noticing, because it means there are two distinct failure modes, and robust code has to anticipate both: the solver exhausts its budget, caught by the converged flag of §3; or the arithmetic overflows first, caught only by an exception handler. Guarding against one and not the other still leaves you exposed, which is why the except clause above names a tuple of exception types rather than RuntimeError alone.

Either way the function refuses to return a number, which is correct. Had irr returned root unconditionally, and had the overflow not intervened, it would have handed back whatever iterate the loop stopped on: a finite, plausible-looking rate for a project with none. An IRR exists only when the discounted-value function actually crosses zero, a condition Descartes' rule of signs supplies in fd01, Exercise 3. A solver cannot manufacture a root that does not exist; what it can do, if you let it, is return a number anyway.

Solution to Exercise 2: Configuration forwarding with **kwargs

def compare_solvers(f, df, a, b, **kwargs):
    """Solve f = 0 by Newton and by bisection, forwarding kwargs to both."""
    x_newton, k_newton, ok_newton = mec_numerical.newton(f, df, 0.5 * (a + b), **kwargs)
    x_bisect, k_bisect, ok_bisect = mec_numerical.bisect(f, a, b, **kwargs)
    if not (ok_newton and ok_bisect):
        raise RuntimeError("at least one solver failed to converge")
    return {"newton": x_newton, "newton_iters": k_newton,
            "bisect": x_bisect, "bisect_iters": k_bisect,
            "gap": abs(x_newton - x_bisect)}


sqrt2 = compare_solvers(lambda x: x ** 2 - 2, lambda x: 2 * x, 1.0, 2.0, tol=1e-12)
print("f(x) = x^2 - 2 on [1, 2]")
for key, value in sqrt2.items():
    print(f"  {key:<14} {value}")

f(x) = x^2 - 2 on [1, 2]
  newton         1.4142135623730951
  newton_iters   4
  bisect         1.4142135623733338
  bisect_iters   37
  gap            2.3869795029440866e-13

solow = compare_solvers(lambda k: s * k ** alpha - delta * k,
                        lambda k: s * alpha * k ** (alpha - 1) - delta,
                        1.0, 20.0, tol=1e-12)
print("Solow: g(k) = s k^alpha - delta k on [1, 20]")
for key, value in solow.items():
    print(f"  {key:<14} {value}")

tol = 1e-9
print(f"\nboth gaps below {tol:.0e}:", sqrt2["gap"] < tol and solow["gap"] < tol)
assert sqrt2["gap"] < tol and solow["gap"] < tol
print(f"and both agree with the closed-form k* = {k_star_closed:.12f}: "
      f"{abs(solow['newton'] - k_star_closed) < tol}")

Solow: g(k) = s k^alpha - delta k on [1, 20]
  newton         7.245789314111252
  newton_iters   4
  bisect         7.245789314103604
  bisect_iters   36
  gap            7.648104372037778e-12

both gaps below 1e-09: True
and both agree with the closed-form k* = 7.245789314111: True

Which takes more iterations, and does that settle it? Bisection takes far more, around forty, against a handful for Newton, because it halves the bracket each step, gaining a fixed number of bits per iteration, while Newton doubles the number of correct digits.

But iteration counts are not the currency that matters, which is the point of Exercise 3. Each bisection step costs one evaluation of $f$; each Newton step costs one of $f$ and one of $f'$. When $f$ is expensive, a nested fixed point, a simulation, a likelihood, the comparison must be made in evaluations, and Newton's advantage narrows. It does not usually reverse on a smooth problem like this one, but the reason to prefer Newton is its convergence rate, not its iteration count, and the reason to keep bisection is that it is guaranteed to converge on any bracket with a sign change, whereas Newton can diverge from a bad start. The two solvers are complements: bisection to get near, Newton to finish.

Solution to Exercise 3: A closure that counts

def counted(func):
    """Wrap func to count its calls. Returns (wrapped, get_count)."""
    calls = 0

    def wrapped(x):
        nonlocal calls
        calls += 1
        return func(x)

    return wrapped, lambda: calls


f_raw, df_raw = lambda x: x ** 2 - 2, lambda x: 2 * x

f_n, count_f_n = counted(f_raw)
df_n, count_df_n = counted(df_raw)
root_n, iters_n, _ = mec_numerical.newton(f_n, df_n, 1.5, tol=1e-12)

f_b, count_f_b = counted(f_raw)
root_b, iters_b, _ = mec_numerical.bisect(f_b, 1.0, 2.0, tol=1e-12)

print(f"{'method':<12}{'iterations':>12}{'f evals':>10}{'df evals':>10}{'total':>8}")
print("-" * 52)
print(f"{'Newton':<12}{iters_n:>12}{count_f_n():>10}{count_df_n():>10}"
      f"{count_f_n() + count_df_n():>8}")
print(f"{'bisection':<12}{iters_b:>12}{count_f_b():>10}{0:>10}{count_f_b():>8}")

print(f"\nboth roots agree: {abs(root_n - root_b):.2e} (tolerance 1e-9)")
assert abs(root_n - root_b) < 1e-9

method        iterations   f evals  df evals   total
----------------------------------------------------
Newton                 4         5         4       9
bisection             37        40         0      40

both roots agree: 2.39e-13 (tolerance 1e-9)

Which is cheaper? Newton, but by less than the iteration counts advertise. Measured in iterations the margin is 4 against 37, nearly ten to one; measured honestly in evaluations it is 9 against 40, a little over four to one, because each Newton step buys its speed by paying for two evaluations instead of one. The advantage remains decisive, since quadratic convergence needs so few steps that paying double for each is a bargain. But the factor is not the one the iteration counter reports, and that is the general argument for measuring cost in the units that actually cost something: the subject of fd04.

If $f'$ were unavailable, the arithmetic changes. Approximating it by a one-sided finite difference costs an extra evaluation of $f$ per step, so a Newton step costs two $f$-evaluations rather than one $f$ and one $f'$: cheaper than it looks if $f'$ is expensive, but now carrying the step-size problem measured in fd01, Exercise 4: the derivative is accurate to about half the available digits, which degrades the convergence from quadratic to roughly superlinear. The standard answer is the secant method, which reuses the previous evaluation and costs one $f$ per step with order $\approx 1.618$; fd09 reaches for Brent's method, which combines bisection's guarantee with superlinear speed. In fd10 the derivative becomes available exactly and at negligible cost by automatic differentiation, at which point the trade-off examined here largely disappears.

Solution to Exercise 4: The Solow convergence rate

Claim. For $T(k) = s k^\alpha + (1-\delta)k$ with $0 < \alpha < 1$, $s, \delta > 0$, and $k^\ast$ the positive steady state,

$$ T'(k^\ast) = 1 - \delta(1-\alpha). $$

Proof. Differentiating, $T'(k) = s\alpha k^{\alpha-1} + 1 - \delta$. The steady-state condition $T(k^\ast) = k^\ast$ reduces to

$$ s (k^\ast)^{\alpha} = \delta k^\ast, $$

and dividing both sides by $k^\ast > 0$ gives $s (k^\ast)^{\alpha - 1} = \delta$. Substituting this into the derivative,

$$ T'(k^\ast) = \alpha \cdot s (k^\ast)^{\alpha-1} + 1 - \delta = \alpha\delta + 1 - \delta = 1 - \delta(1-\alpha). \qquad \blacksquare $$

The steady-state condition is used exactly once, to replace $s(k^\ast)^{\alpha-1}$ by $\delta$, which is what removes $s$ from the expression.

Half-life. Since $k_t - k^\ast \approx \lambda^t (k_0 - k^\ast)$ with $\lambda = T'(k^\ast) \in (0,1)$, the error halves when $\lambda^t = 1/2$, that is after

$$ t_{1/2} = \frac{\ln 2}{-\ln\lambda} = \frac{\ln 2}{-\ln\left(1 - \delta(1-\alpha)\right)} \text{ periods.} $$

Why $s$ does not appear. The saving rate enters $T$ only through the level of the production term, and the steady-state condition ties that level to $\delta$: raising $s$ raises $k^\ast$ until the marginal product falls enough that depreciation again absorbs all investment. At whatever steady state results, $s(k^\ast)^{\alpha-1}$ equals $\delta$ by construction. So $s$ relocates the fixed point but leaves the slope of $T$ there unchanged: it moves the destination, not the speed of the journey. Only the curvature of the production function, $\alpha$, and the rate at which capital is lost, $\delta$, govern the speed.

import math

print(f"{'alpha':>7}{'delta':>8}{'predicted':>13}{'measured':>13}{'gap':>10}"
      f"{'half-life':>12}{'periods':>10}")
print("-" * 73)
worst = 0.0
for alpha_ in (0.3, 0.5, 0.7):
    for delta_ in (0.02, 0.05, 0.10):
        step = make_solow_step(s, alpha_, delta_)
        k_star = (s / delta_) ** (1 / (1 - alpha_))
        lam_hat = 1 - delta_ * (1 - alpha_)

        # Iterate until the RELATIVE error is small enough that the linearization
        # holds, but not so small that cancellation destroys the ratio. A fixed
        # iteration budget will not do: the slow cases need thousands of periods.
        k, periods = 0.1, 0
        while abs(k - k_star) / k_star > 1e-7 and periods < 20_000:
            k = step(k)
            periods += 1
        assert periods < 20_000, f"({alpha_}, {delta_}) never reached the asymptotic regime"

        measured = (step(k) - k_star) / (k - k_star)
        gap = abs(measured - lam_hat)
        worst = max(worst, gap)
        half_life = math.log(2) / (-math.log(lam_hat))
        print(f"{alpha_:>7.1f}{delta_:>8.2f}{lam_hat:>13.6f}{measured:>13.6f}"
              f"{gap:>10.1e}{half_life:>12.1f}{periods:>10}")

tol = 1e-6
print("-" * 73)
print(f"worst gap across the nine cases = {worst:.1e}   (tolerance {tol:.0e})")
assert worst < tol, "the measured rates do not match the prediction"
print("check passed.")

  alpha   delta    predicted     measured       gap   half-life   periods
-------------------------------------------------------------------------
    0.3    0.02     0.986000     0.986000   1.1e-10        49.2      1168
    0.3    0.05     0.965000     0.965000   5.0e-10        19.5       462
    0.3    0.10     0.930000     0.930000   1.8e-11         9.6       227
    0.5    0.02     0.990000     0.990000   1.0e-09        69.0      1672
    0.5    0.05     0.975000     0.975000   1.1e-09        27.4       663
    0.5    0.10     0.950000     0.950000   1.4e-09        13.5       326
    0.7    0.02     0.994000     0.994000   8.6e-10       115.2      2874
    0.7    0.05     0.985000     0.985000   7.8e-10        45.9      1140
    0.7    0.10     0.970000     0.970000   9.2e-10        22.8       561
-------------------------------------------------------------------------
worst gap across the nine cases = 1.4e-09   (tolerance 1e-06)
check passed.

The measured ratios track the prediction across all nine parameter pairs. The tolerance here is looser than in §9, $10^{-4}$ rather than $10^{-6}$, and honestly so: the loop stops as soon as the error falls below $10^{-8}$, which is not always deep enough into the asymptotic regime for the linearization to hold to six decimals. Choosing a tolerance means being explicit about how good an approximation you have actually earned, and quoting six digits here would be quoting the linearization's error as though it were the code's.

Which case looks like the data? The empirical convergence rate of roughly $2\%$ per year corresponds to $\lambda \approx 0.98$ and a half-life near 35 years. No row hits it exactly, but two bracket it at a conventional depreciation rate: $\alpha = 0.5$, $\delta = 0.05$ gives a half-life of 27 years, and $\alpha = 0.7$, $\delta = 0.05$ gives 46. Solving (9.3) directly, $\delta(1-\alpha) = 0.02$ at $\delta = 0.05$ requires

$$ \alpha = 1 - \frac{0.02}{0.05} = 0.6 . $$

So matching the observed speed of convergence demands a capital share near $0.6$, twice the $0.3$ of the national accounts. That is the classic tension, and its classic resolution is to read $k$ as a broad aggregate including human capital: precisely the argument of Mankiw, Romer and Weil (1992). The alternative route, holding $\alpha = 0.3$ and lowering $\delta$ to $0.029$, is not credible as a depreciation rate.

The methodological point is the one worth keeping. The model does not merely fail to match the data at textbook parameters; through (9.3) it states exactly which parameter must move and by how much. A model that quantifies its own discrepancy is more useful than one that merely fits.

Solution to Exercise 5: A @dataclass for a portfolio

from dataclasses import dataclass, field

@dataclass
class Position:
    ticker: str
    quantity: int
    price: float

    def value(self):
        return self.quantity * self.price


@dataclass
class Portfolio:
    positions: list = field(default_factory=list)     # NOT positions: list = []

    def total_value(self):
        return sum(p.value() for p in self.positions)

    def weights(self):
        total = self.total_value()
        return {p.ticker: p.value() / total for p in self.positions}


portfolio = Portfolio([
    Position("AAPL", 100, 189.50),
    Position("MSFT", 50, 412.30),
    Position("BRK.B", 20, 405.75),
])

print(portfolio, "\n")
for ticker, weight in portfolio.weights().items():
    print(f"  {ticker:<7} {weight:>7.2%}")

Portfolio(positions=[Position(ticker='AAPL', quantity=100, price=189.5), Position(ticker='MSFT', quantity=50, price=412.3), Position(ticker='BRK.B', quantity=20, price=405.75)]) 

  AAPL     39.74%
  MSFT     43.24%
  BRK.B    17.02%

# route 1: sum the method over positions;  route 2: a comprehension over raw fields
total_method = portfolio.total_value()
total_raw = sum(p.quantity * p.price for p in portfolio.positions)

gap, tol = abs(total_method - total_raw), 1e-10
print(f"total via value()      = {total_method:>12,.2f}")
print(f"total via raw fields   = {total_raw:>12,.2f}")
print(f"gap = {gap:.1e}   (tolerance {tol:.0e})")
assert gap < tol

weight_sum = sum(portfolio.weights().values())
print(f"\nweights sum to {weight_sum:.16f}, |sum - 1| = {abs(weight_sum - 1):.1e}"
      f"   (tolerance 1e-12)")
assert abs(weight_sum - 1) < 1e-12
print("check passed.")

total via value()      =    47,680.00
total via raw fields   =    47,680.00
gap = 0.0e+00   (tolerance 1e-10)

weights sum to 1.0000000000000000, |sum - 1| = 0.0e+00   (tolerance 1e-12)
check passed.

The design question. If two portfolios were constructed from the same list object, they would share it, fd02 §7, and appending a position to one would silently add it to the other. Both would then report the same, wrong, total.

Note that the dataclass field is written positions: list = field(default_factory=list) rather than positions: list = []. This is the mutable-default trap of §3 in its dataclass form, and it is one of the few places where Python protects you: a mutable default in a dataclass raises ValueError at class-definition time rather than failing silently at runtime. default_factory calls list() afresh for each instance, which is the dataclass spelling of the None-sentinel idiom.

Two remarks on the weights check. First, it is nearly free but not vacuous: it would catch a total computed over a different set of positions than the one iterated, which is exactly the failure the aliasing bug produces. Second, the sum of the three weights is not exactly $1$ in binary floating point, hence the tolerance, for precisely the reason fd01 §5 gave for never comparing floats with ==. A shares-of-a-total check is one of the most common places where that rule gets forgotten.