Tutorial: Keen Predator-Prey on the Pacioli Manifold¶

Steve Keen's 1995 model is the canonical example of Minsky dynamics: the inherent instability of a capitalist economy driven by private debt accumulation. It is a three-variable dynamical system:

Variable	Symbol	Meaning
Wage share	ω	Labour's share of GDP
Employment rate	λ	Fraction of workforce employed
Debt ratio	d	Private debt / GDP

The dynamics are:

$$\dot{\omega} = \omega\,(\Phi(\lambda) - \alpha)$$ $$\dot{\lambda} = \lambda\,(\kappa(\pi)/\nu - \alpha - \beta - \gamma)$$ $$\dot{d} = \kappa(\pi) - \pi - (\alpha + \beta)\,d$$

where $\pi = 1 - \omega - r d$ is the profit share, $\Phi(\lambda)$ is the Phillips curve, and $\kappa(\pi)$ is the investment function.

What this tutorial demonstrates:

1. Simulate the Keen ODE using econiac.economics.minsky
2. Inspect sectoral balance sheets on the Pacioli manifold
3. Model the Minsky moment as a PCL circuit: choose(β, invest, deleverage)
4. Route between strategies using TIR thermodynamic routing
5. Attribute the crisis using thermal Shapley values
6. Calibrate against noisy national accounts using conservation_loss

---

Reference: Keen (1995) Finance and Economic Breakdown.
EconIAC implementation: Buckley (2026), doi:10.5281/zenodo.20262070

In [2]:

# !pip install econiac "jax[cpu]"

# !pip install econiac "jax[cpu]"

In [3]:

# Install econiac if running in Colab
# !pip install econiac jax[cpu]

import jax.numpy as jnp
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

from econiac.core.manifold import BalanceSheet, add_residual_sector, residual_magnitude, three_sector_sfc
from econiac.economics.minsky import keen_simulate
from econiac.pcl import flow, sequence, choose, conservation_loss, compile as pcl_compile
from econiac.pcl.combinators import Computation
from econiac.routing import tir_from_scores, route
from econiac.routing.attribution import full_shapley_analysis

plt.rcParams.update({"figure.dpi": 120, "font.size": 11})

# Install econiac if running in Colab # !pip install econiac jax[cpu] import jax.numpy as jnp import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from econiac.core.manifold import BalanceSheet, add_residual_sector, residual_magnitude, three_sector_sfc from econiac.economics.minsky import keen_simulate from econiac.pcl import flow, sequence, choose, conservation_loss, compile as pcl_compile from econiac.pcl.combinators import Computation from econiac.routing import tir_from_scores, route from econiac.routing.attribution import full_shapley_analysis plt.rcParams.update({"figure.dpi": 120, "font.size": 11})

1. Simulate the Keen ODE¶

keen_simulate integrates the three-variable system by Euler steps. The state vector is $[\omega, \lambda, d]$ at each timestep.

In [4]:

times, traj = keen_simulate(omega0=0.80, lam0=0.94, d0=0.10, T=150.0, dt=0.05)

omega = traj[:, 0]
lam   = traj[:, 1]
d     = traj[:, 2]

print(f"Steps: {len(times)}  |  T = {float(times[-1]):.0f} years")
print(f"Initial:  ω={float(omega[0]):.3f}  λ={float(lam[0]):.3f}  d={float(d[0]):.3f}")
print(f"Final:    ω={float(omega[-1]):.3f}  λ={float(lam[-1]):.3f}  d={float(d[-1]):.3f}")
print(f"Peak debt: {float(d.max()):.3f}  at t = {float(times[int(d.argmax())]):.1f}")

times, traj = keen_simulate(omega0=0.80, lam0=0.94, d0=0.10, T=150.0, dt=0.05) omega = traj[:, 0] lam = traj[:, 1] d = traj[:, 2] print(f"Steps: {len(times)} | T = {float(times[-1]):.0f} years") print(f"Initial: ω={float(omega[0]):.3f} λ={float(lam[0]):.3f} d={float(d[0]):.3f}") print(f"Final: ω={float(omega[-1]):.3f} λ={float(lam[-1]):.3f} d={float(d[-1]):.3f}") print(f"Peak debt: {float(d.max()):.3f} at t = {float(times[int(d.argmax())]):.1f}")

Steps: 3001  |  T = 150 years
Initial:  ω=0.800  λ=0.940  d=0.100
Final:    ω=0.060  λ=0.030  d=81.884
Peak debt: 93.304  at t = 116.2

In [5]:

fig, axes = plt.subplots(1, 3, figsize=(13, 4))
t = times

axes[0].plot(t, omega, color="steelblue")
axes[0].set_xlabel("Years"); axes[0].set_ylabel("ω"); axes[0].set_title("Wage share ω(t)")
axes[0].axhline(float(omega[0]), ls="--", color="grey", lw=0.8, label="initial")
axes[0].legend(fontsize=9)

axes[1].plot(t, lam, color="seagreen")
axes[1].set_xlabel("Years"); axes[1].set_ylabel("λ"); axes[1].set_title("Employment rate λ(t)")

axes[2].plot(t, d, color="crimson")
axes[2].axvline(float(times[int(d.argmax())]), ls="--", color="grey", lw=0.8, label="peak debt")
axes[2].set_xlabel("Years"); axes[2].set_ylabel("d"); axes[2].set_title("Debt ratio d(t)")
axes[2].legend(fontsize=9)

fig.suptitle("Keen predator-prey dynamics", fontsize=13, y=1.02)
plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 3, figsize=(13, 4)) t = times axes[0].plot(t, omega, color="steelblue") axes[0].set_xlabel("Years"); axes[0].set_ylabel("ω"); axes[0].set_title("Wage share ω(t)") axes[0].axhline(float(omega[0]), ls="--", color="grey", lw=0.8, label="initial") axes[0].legend(fontsize=9) axes[1].plot(t, lam, color="seagreen") axes[1].set_xlabel("Years"); axes[1].set_ylabel("λ"); axes[1].set_title("Employment rate λ(t)") axes[2].plot(t, d, color="crimson") axes[2].axvline(float(times[int(d.argmax())]), ls="--", color="grey", lw=0.8, label="peak debt") axes[2].set_xlabel("Years"); axes[2].set_ylabel("d"); axes[2].set_title("Debt ratio d(t)") axes[2].legend(fontsize=9) fig.suptitle("Keen predator-prey dynamics", fontsize=13, y=1.02) plt.tight_layout() plt.show()

In [6]:

# Phase portrait: (λ, ω) — the classic Keen portrait
fig, ax = plt.subplots(figsize=(6, 5))
sc = ax.scatter(lam, omega, c=t, cmap="plasma", s=2)
ax.set_xlabel("Employment rate λ"); ax.set_ylabel("Wage share ω")
ax.set_title("Phase portrait (λ, ω)\ncolour = time")
plt.colorbar(sc, ax=ax, label="Years")
plt.tight_layout()
plt.show()

# Phase portrait: (λ, ω) — the classic Keen portrait fig, ax = plt.subplots(figsize=(6, 5)) sc = ax.scatter(lam, omega, c=t, cmap="plasma", s=2) ax.set_xlabel("Employment rate λ"); ax.set_ylabel("Wage share ω") ax.set_title("Phase portrait (λ, ω)\ncolour = time") plt.colorbar(sc, ax=ax, label="Years") plt.tight_layout() plt.show()

2. Balance sheets on the Pacioli manifold¶

Every point on the Keen trajectory corresponds to a sectoral balance sheet. The conservation law $\partial^2 = 0$ — every asset is someone else's liability — must hold at all times.

We can also simulate measurement error: real national accounts always have a statistical discrepancy. add_residual_sector absorbs any imbalance into a designated _residual sector, making the extended balance sheet exactly flat while keeping the discrepancy visible and measurable.

In [7]:

sectors     = ["workers", "firms", "banks"]
instruments = ["wage_share", "debt_ratio"]

def make_bs(i):
    w, dv = float(omega[i]), float(d[i])
    return BalanceSheet(
        positions=jnp.array([[w, 0.0], [-w, -dv], [0.0, dv]]),
        sectors=sectors, instruments=instruments,
    )

for label, idx in [("t=0 (initial)", 0),
                   (f"t={float(times[int(d.argmax())]):.0f} (peak debt)", int(d.argmax())),
                   ("t=150 (final)", -1)]:
    bs = make_bs(idx)
    print(f"{label:25s}  {bs}  col_sums={bs.column_sums()}")

sectors = ["workers", "firms", "banks"] instruments = ["wage_share", "debt_ratio"] def make_bs(i): w, dv = float(omega[i]), float(d[i]) return BalanceSheet( positions=jnp.array([[w, 0.0], [-w, -dv], [0.0, dv]]), sectors=sectors, instruments=instruments, ) for label, idx in [("t=0 (initial)", 0), (f"t={float(times[int(d.argmax())]):.0f} (peak debt)", int(d.argmax())), ("t=150 (final)", -1)]: bs = make_bs(idx) print(f"{label:25s} {bs} col_sums={bs.column_sums()}")

t=0 (initial)              BalanceSheet(3 sectors × 2 instruments, ✓ consistent)  col_sums=[0. 0.]
t=116 (peak debt)          BalanceSheet(3 sectors × 2 instruments, ✓ consistent)  col_sums=[0. 0.]
t=150 (final)              BalanceSheet(3 sectors × 2 instruments, ✓ consistent)  col_sums=[0. 0.]

In [8]:

# Track the residual magnitude (= statistical discrepancy) along the trajectory
# Simulate a 3% survey measurement error in the wage share
residuals = []
for i in range(0, len(times), 10):
    w, dv = float(omega[i]), float(d[i])
    noisy = jnp.array([[w + 0.03, 0.0], [-w, -dv], [0.0, dv]])
    bs_noisy = BalanceSheet(positions=noisy, sectors=sectors, instruments=instruments)
    bs_rec   = add_residual_sector(bs_noisy)
    residuals.append(residual_magnitude(bs_rec))

fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(times[::10], residuals, color="darkorange")
ax.set_xlabel("Years"); ax.set_ylabel("||residual||")
ax.set_title("Statistical discrepancy magnitude along the Keen trajectory\n(3% survey error in wage share)")
plt.tight_layout()
plt.show()

# Track the residual magnitude (= statistical discrepancy) along the trajectory # Simulate a 3% survey measurement error in the wage share residuals = [] for i in range(0, len(times), 10): w, dv = float(omega[i]), float(d[i]) noisy = jnp.array([[w + 0.03, 0.0], [-w, -dv], [0.0, dv]]) bs_noisy = BalanceSheet(positions=noisy, sectors=sectors, instruments=instruments) bs_rec = add_residual_sector(bs_noisy) residuals.append(residual_magnitude(bs_rec)) fig, ax = plt.subplots(figsize=(8, 3)) ax.plot(times[::10], residuals, color="darkorange") ax.set_xlabel("Years"); ax.set_ylabel("||residual||") ax.set_title("Statistical discrepancy magnitude along the Keen trajectory\n(3% survey error in wage share)") plt.tight_layout() plt.show()

3. The Minsky moment as a PCL circuit¶

The core of the Keen model is the investment switch: when profit falls, firms must choose between borrowing more (Ponzi finance) and paying down debt (deleveraging).

In PCL, this is the transistor:

choose(β, invest, deleverage)

• β = 0: equal weight to both strategies — the system hedges
• β → ∞: all weight on whichever strategy produces higher net worth — a hard switch

This is not a metaphor. At β → ∞, choose recovers the tropical (max, +) semiring — the mathematical definition of digital switching. At finite β it is the Gibbs ensemble average: a smooth, differentiable interpolation between exploration and exploitation.

In [9]:

bs_pcl = BalanceSheet(
    positions=jnp.array([[100.0, 0.0], [0.0, -80.0], [-100.0, 80.0]]),
    sectors=["households", "firms", "banks"],
    instruments=["deposits", "loans"],
)

invest     = flow("banks", "firms",      "loans", 20.0)
deleverage = flow("firms", "banks",      "loans", 20.0)
wages      = flow("firms", "households", "deposits", 30.0)
taxes      = flow("households", "firms", "deposits",  5.0)

# Sweep β from 0 (pure hedge) to 10 (near-digital switch)
betas       = [0.0, 0.5, 1.0, 2.0, 5.0, 10.0]
firm_loans  = []

for beta in betas:
    circuit = sequence(wages, sequence(taxes, choose(beta, invest, deleverage)))
    result  = circuit(bs_pcl)
    firm_loans.append(float(result.positions[1, 1]))

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(betas, firm_loans, "o-", color="steelblue", lw=2)
ax.axhline(-60, ls="--", color="crimson",  lw=1, label="full invest (loans = −60)")
ax.axhline(-100, ls="--", color="seagreen", lw=1, label="full deleverage (loans = −100)")
ax.set_xlabel("β (rationality / commitment)")
ax.set_ylabel("Firm loans after quarterly circuit")
ax.set_title("Minsky moment: choose(β, invest, deleverage)\nβ→∞ commits to higher-value strategy")
ax.legend()
plt.tight_layout()
plt.show()

bs_pcl = BalanceSheet( positions=jnp.array([[100.0, 0.0], [0.0, -80.0], [-100.0, 80.0]]), sectors=["households", "firms", "banks"], instruments=["deposits", "loans"], ) invest = flow("banks", "firms", "loans", 20.0) deleverage = flow("firms", "banks", "loans", 20.0) wages = flow("firms", "households", "deposits", 30.0) taxes = flow("households", "firms", "deposits", 5.0) # Sweep β from 0 (pure hedge) to 10 (near-digital switch) betas = [0.0, 0.5, 1.0, 2.0, 5.0, 10.0] firm_loans = [] for beta in betas: circuit = sequence(wages, sequence(taxes, choose(beta, invest, deleverage))) result = circuit(bs_pcl) firm_loans.append(float(result.positions[1, 1])) fig, ax = plt.subplots(figsize=(7, 4)) ax.plot(betas, firm_loans, "o-", color="steelblue", lw=2) ax.axhline(-60, ls="--", color="crimson", lw=1, label="full invest (loans = −60)") ax.axhline(-100, ls="--", color="seagreen", lw=1, label="full deleverage (loans = −100)") ax.set_xlabel("β (rationality / commitment)") ax.set_ylabel("Firm loans after quarterly circuit") ax.set_title("Minsky moment: choose(β, invest, deleverage)\nβ→∞ commits to higher-value strategy") ax.legend() plt.tight_layout() plt.show()

4. TIR routing at the Minsky moment¶

At each moment the economy faces three options: invest (expand), deleverage (contract), or hold steady. The Thermodynamic Information Routing (TIR) framework weights these options by a Gibbs distribution over their expected profit scores.

• Low β: near-uniform weights — the system explores all options
• High β: concentrated weights — the system exploits the best option

In [10]:

labels        = ["invest", "deleverage", "hold"]
profit_scores = [0.15, 0.25, 0.20]   # deleverage has highest expected return at Minsky moment

beta_sweep = [0.2, 0.5, 1.0, 2.0, 5.0, 10.0]
all_weights = []
for b in beta_sweep:
    w = route(tir_from_scores(labels, profit_scores, beta=b))
    all_weights.append([float(wi) for wi in w])

all_weights = list(zip(*all_weights))   # transpose: (n_options, n_betas)

fig, ax = plt.subplots(figsize=(8, 4))
colors = ["crimson", "seagreen", "steelblue"]
for i, (label, wts, c) in enumerate(zip(labels, all_weights, colors)):
    ax.plot(beta_sweep, wts, "o-", color=c, label=label, lw=2)

ax.set_xlabel("β"); ax.set_ylabel("Routing weight")
ax.set_title("TIR routing at the Minsky moment\nDeleverage wins as β increases")
ax.legend(); ax.set_ylim(0, 1)
plt.tight_layout()
plt.show()

labels = ["invest", "deleverage", "hold"] profit_scores = [0.15, 0.25, 0.20] # deleverage has highest expected return at Minsky moment beta_sweep = [0.2, 0.5, 1.0, 2.0, 5.0, 10.0] all_weights = [] for b in beta_sweep: w = route(tir_from_scores(labels, profit_scores, beta=b)) all_weights.append([float(wi) for wi in w]) all_weights = list(zip(*all_weights)) # transpose: (n_options, n_betas) fig, ax = plt.subplots(figsize=(8, 4)) colors = ["crimson", "seagreen", "steelblue"] for i, (label, wts, c) in enumerate(zip(labels, all_weights, colors)): ax.plot(beta_sweep, wts, "o-", color=c, label=label, lw=2) ax.set_xlabel("β"); ax.set_ylabel("Routing weight") ax.set_title("TIR routing at the Minsky moment\nDeleverage wins as β increases") ax.legend(); ax.set_ylim(0, 1) plt.tight_layout() plt.show()

5. Thermal Shapley attribution¶

Which sector is most responsible for the crisis — and most capable of averting it?

Thermal Shapley values generalise the classical Shapley value by weighting coalitions according to the Gibbs distribution. The β parameter controls how much weight is placed on the highest-value coalitions:

• β = 0: classical Shapley values (uniform weights)
• β → ∞: all weight on the best coalition (bottleneck player = single point of leverage)

The value function here measures how much each coalition could improve GDP by acting prudently (moderating wages, deleveraging, tightening credit).

In [11]:

def crisis_value(coalition: frozenset) -> float:
    base = -float(d.max())
    if 0 in coalition: base += 0.15   # workers moderate wages
    if 1 in coalition: base += 0.35   # firms deleverage
    if 2 in coalition: base += 0.20   # banks tighten credit
    return base

player_names = ["workers", "firms", "banks"]
beta_sweep_shapley = [0.0, 0.5, 1.0, 2.0, 5.0]
all_phi = []

for b in beta_sweep_shapley:
    result, _ = full_shapley_analysis(crisis_value, n_players=3, beta=b)
    all_phi.append([float(v) for v in result.phi])

all_phi = list(zip(*all_phi))

fig, ax = plt.subplots(figsize=(8, 4))
colors = ["steelblue", "crimson", "seagreen"]
for name, phis, c in zip(player_names, all_phi, colors):
    ax.plot(beta_sweep_shapley, phis, "o-", color=c, label=name, lw=2)

ax.set_xlabel("β"); ax.set_ylabel("Shapley value φ")
ax.set_title("Thermal Shapley values by sector\nFirms are the bottleneck — largest leverage for crisis mitigation")
ax.legend()
plt.tight_layout()
plt.show()

# Bar chart at β = 2
result_b2, _ = full_shapley_analysis(crisis_value, n_players=3, beta=2.0)
fig, ax = plt.subplots(figsize=(5, 3.5))
bars = ax.bar(player_names, [float(v) for v in result_b2.phi], color=colors)
ax.set_ylabel("Shapley value φ (β=2)")
ax.set_title(f"Crisis attribution at β=2\nBottleneck: {player_names[result_b2.bottleneck_player]}")
plt.tight_layout()
plt.show()

def crisis_value(coalition: frozenset) -> float: base = -float(d.max()) if 0 in coalition: base += 0.15 # workers moderate wages if 1 in coalition: base += 0.35 # firms deleverage if 2 in coalition: base += 0.20 # banks tighten credit return base player_names = ["workers", "firms", "banks"] beta_sweep_shapley = [0.0, 0.5, 1.0, 2.0, 5.0] all_phi = [] for b in beta_sweep_shapley: result, _ = full_shapley_analysis(crisis_value, n_players=3, beta=b) all_phi.append([float(v) for v in result.phi]) all_phi = list(zip(*all_phi)) fig, ax = plt.subplots(figsize=(8, 4)) colors = ["steelblue", "crimson", "seagreen"] for name, phis, c in zip(player_names, all_phi, colors): ax.plot(beta_sweep_shapley, phis, "o-", color=c, label=name, lw=2) ax.set_xlabel("β"); ax.set_ylabel("Shapley value φ") ax.set_title("Thermal Shapley values by sector\nFirms are the bottleneck — largest leverage for crisis mitigation") ax.legend() plt.tight_layout() plt.show() # Bar chart at β = 2 result_b2, _ = full_shapley_analysis(crisis_value, n_players=3, beta=2.0) fig, ax = plt.subplots(figsize=(5, 3.5)) bars = ax.bar(player_names, [float(v) for v in result_b2.phi], color=colors) ax.set_ylabel("Shapley value φ (β=2)") ax.set_title(f"Crisis attribution at β=2\nBottleneck: {player_names[result_b2.bottleneck_player]}") plt.tight_layout() plt.show()

6. Calibration: conservation_loss against noisy data¶

When calibrating a model against real national accounts, we cannot demand that the data satisfy $\partial^2 = 0$ to machine precision — the ONS and BEA publish data with known measurement uncertainty.

conservation_loss(comp, bs, sigma) replaces the binary typecheck with a differentiable penalty:

$$\mathcal{L} = \frac{\|\text{col\_sums}(\text{comp}(bs))\|^2}{\sigma^2}$$

• At σ → 0: infinite penalty (hard conservation)
• At finite σ: Bayesian likelihood under Gaussian measurement noise

The loss is fully differentiable via JAX — gradients flow through the entire circuit including every choose and fold gate.

In [12]:

circuit = sequence(wages, sequence(taxes, choose(2.0, invest, deleverage)))

# Well-specified conserving circuit: loss = 0 regardless of σ
sigmas = [0.1, 0.5, 1.0, 5.0, 10.0, 50.0]
losses_good = [float(conservation_loss(circuit, bs_pcl, sigma=s)) for s in sigmas]

# Non-conserving computation (simulates model misspecification)
def leaky_fn(bs):
    return BalanceSheet(
        positions=bs.positions.at[0, 0].add(10.0),
        sectors=bs.sectors, instruments=bs.instruments,
    )
leaky = Computation(name="leaky", fn=leaky_fn)
losses_bad = [float(conservation_loss(leaky, bs_pcl, sigma=s)) for s in sigmas]

fig, ax = plt.subplots(figsize=(8, 4))
ax.loglog(sigmas, losses_good, "o-", color="seagreen",  label="conserving circuit (loss = 0)")
ax.loglog(sigmas, losses_bad,  "o-", color="crimson",   label="leaky model (10-unit error)")
ax.set_xlabel("σ (measurement noise scale)"); ax.set_ylabel("conservation_loss")
ax.set_title("conservation_loss vs. σ\nConserving circuits: zero loss. Leaky models: penalised by 1/σ²")
ax.legend()
plt.tight_layout()
plt.show()

print("conservation_loss for the well-specified circuit:")
for s, l in zip(sigmas, losses_good):
    print(f"  σ={s:5.1f}  loss={l:.6f}")

circuit = sequence(wages, sequence(taxes, choose(2.0, invest, deleverage))) # Well-specified conserving circuit: loss = 0 regardless of σ sigmas = [0.1, 0.5, 1.0, 5.0, 10.0, 50.0] losses_good = [float(conservation_loss(circuit, bs_pcl, sigma=s)) for s in sigmas] # Non-conserving computation (simulates model misspecification) def leaky_fn(bs): return BalanceSheet( positions=bs.positions.at[0, 0].add(10.0), sectors=bs.sectors, instruments=bs.instruments, ) leaky = Computation(name="leaky", fn=leaky_fn) losses_bad = [float(conservation_loss(leaky, bs_pcl, sigma=s)) for s in sigmas] fig, ax = plt.subplots(figsize=(8, 4)) ax.loglog(sigmas, losses_good, "o-", color="seagreen", label="conserving circuit (loss = 0)") ax.loglog(sigmas, losses_bad, "o-", color="crimson", label="leaky model (10-unit error)") ax.set_xlabel("σ (measurement noise scale)"); ax.set_ylabel("conservation_loss") ax.set_title("conservation_loss vs. σ\nConserving circuits: zero loss. Leaky models: penalised by 1/σ²") ax.legend() plt.tight_layout() plt.show() print("conservation_loss for the well-specified circuit:") for s, l in zip(sigmas, losses_good): print(f" σ={s:5.1f} loss={l:.6f}")

conservation_loss for the well-specified circuit:
  σ=  0.1  loss=0.000000
  σ=  0.5  loss=0.000000
  σ=  1.0  loss=0.000000
  σ=  5.0  loss=0.000000
  σ= 10.0  loss=0.000000
  σ= 50.0  loss=0.000000

Summary¶

This tutorial demonstrated the full econiac stack on the Keen predator-prey model:

Step	Module	What it showed
ODE simulation	econiac.economics.minsky	150-year Keen trajectory; phase portrait
Balance sheets	econiac.core.manifold	∂²=0 conservation at every timestep; residual sector for noisy data
PCL circuit	econiac.pcl	choose(β, invest, deleverage) as the Minsky transistor; β-sweep from hedge to commit
TIR routing	econiac.routing	Three-way Gibbs routing; deleverage wins at high β
Thermal Shapley	econiac.routing.attribution	Firms are the bottleneck (largest crisis-mitigation leverage)
Calibration	econiac.pcl.conservation_loss	Differentiable penalty; 1/σ² scaling with measurement noise

Next: GEMMES — Keen + climate on the Pacioli manifold

---

EconIAC is open source: pip install econiac
Theory: Buckley (2026), Portfolio G papers, doi:10.5281/zenodo.20259495