Tutorial: GEMMES — Keen + Climate on the Pacioli Manifold¶

GEMMES (General Monetary and Multisectoral Macrodynamics for the Ecological Shift) extends the Keen predator-prey model with a carbon cycle and climate damage function. It is the canonical model for studying the interaction between private debt dynamics and climate risk.

The state vector is now five-dimensional:

Variable	Symbol	Meaning
Wage share	ω	Labour's share of GDP
Employment rate	λ	Fraction of workforce employed
Debt ratio	d	Private debt / GDP
Temperature anomaly	ΔT	Global warming above pre-industrial (°C)
CO₂ concentration	CO₂	Atmospheric carbon (GtC)

The climate variables feed back into the economy via the Nordhaus damage function:

$$D(\Delta T) = 1 - \frac{1}{1 + \pi_1 \Delta T + \pi_2 \Delta T^2}$$

which reduces the effective profit share: $\pi_\text{net} = (1 - \omega - rd)(1 - D(\Delta T))$.

What this tutorial adds beyond the Keen tutorial:

1. Climate-augmented ODE with CO₂ and temperature dynamics
2. CurvedBalanceSheet — stranded asset risk as non-zero curvature field F
3. PCL fold(β, [brown, green, deleverage]) — the three-way green transition switch
4. TIR routing: how a carbon tax shifts investment weights
5. 4-player thermal Shapley (workers / firms / banks / policy)

---

Reference: Bovari, Giraud, McIsaac & Chancel (2018) Ecological Economics 147, 383–398.
EconIAC implementation: Buckley (2026), doi:10.5281/zenodo.20262070

In [1]:

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

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

from econiac.core.manifold import (
    BalanceSheet, CurvedBalanceSheet, holonomy,
)
from econiac.pcl import flow, sequence, fold, conservation_loss, compile as pcl_compile
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 import jax.numpy as jnp import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from econiac.core.manifold import ( BalanceSheet, CurvedBalanceSheet, holonomy, ) from econiac.pcl import flow, sequence, fold, conservation_loss, compile as pcl_compile 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. The GEMMES ODE system¶

We implement the five-variable system directly in JAX so it is differentiable end-to-end. The key additions over the Keen model are:

• Nordhaus damage $D(\Delta T)$ erodes the profit share
• Radiative forcing $F(\text{CO}_2) = F_{2\times} \log_2(\text{CO}_2 / \text{CO}_{2,\text{pre}})$ drives warming
• Carbon cycle $\dot{\text{CO}}_2 = \sigma(1-n)Y - \delta_{\text{CO}_2} \cdot \text{CO}_2$ where $n$ is the abatement fraction

In [2]:

PARAMS = dict(
    alpha=0.02, beta=0.01, gamma=0.01, nu=3.0, r=0.03,
    phi0=-0.04, phi1=0.0006,
    kappa0=0.0, kappa1=0.025, kappa2=2.0,   # damped investment for stable orbit
    pi1=0.0, pi2=0.00236,                    # Nordhaus (2008) damage coefficients
    tau_T=50.0, lambda_T=1.13, F2x=3.7,
    CO2_pre=280.0, sigma=0.02, delta_CO2=0.01, n_tax=0.10,
)

def damage(T, p):
    return 1.0 - 1.0 / (1.0 + p["pi1"] * T + p["pi2"] * T**2)

def gemmes_ode(state, t, p):
    omega, lam, d, T, CO2 = state
    D     = damage(T, p)
    pi    = (1.0 - omega - p["r"] * d) * (1.0 - D)
    phi   = p["phi0"] + p["phi1"] / jnp.maximum(1.0 - lam, 0.01)**2
    kappa = (p["kappa0"] + p["kappa1"] * jnp.exp(p["kappa2"] * pi)) / p["nu"]
    F     = p["F2x"] * jnp.log(jnp.maximum(CO2, 1.0) / p["CO2_pre"]) / jnp.log(2.0)
    return jnp.array([
        omega * (phi - p["alpha"]),
        lam   * (kappa - p["alpha"] - p["beta"]),
        kappa - pi - (p["alpha"] + p["beta"]) * d,
        (F - p["lambda_T"] * T) / p["tau_T"],
        p["sigma"] * (1.0 - p["n_tax"]) * 100.0 - p["delta_CO2"] * CO2,
    ])

def simulate(omega0=0.78, lam0=0.95, d0=0.05, T0=0.85, CO2_0=395.0,
             T_end=80.0, dt=0.02, params=None):
    p = params or PARAMS
    state = jnp.array([omega0, lam0, d0, T0, CO2_0])
    n     = int(round(T_end / dt))
    times, traj = [0.0], [state]
    t = 0.0
    for _ in range(n):
        state = state + gemmes_ode(state, t, p) * dt
        t    += dt
        times.append(t); traj.append(state)
    return jnp.array(times), jnp.stack(traj)

times, traj = simulate()
omega, lam, d, T_anom, CO2 = [traj[:, i] for i in range(5)]

print(f"T = {float(times[-1]):.0f} years  |  {len(times)} steps")
print(f"Peak ΔT:   {float(T_anom.max()):.2f}°C  at t={float(times[int(T_anom.argmax())]):.0f}")
print(f"Peak debt: {float(d.max()):.3f}     at t={float(times[int(d.argmax())]):.0f}")
print(f"Damage at peak ΔT: {float(damage(T_anom.max(), PARAMS))*100:.2f}% of GDP")

PARAMS = dict( alpha=0.02, beta=0.01, gamma=0.01, nu=3.0, r=0.03, phi0=-0.04, phi1=0.0006, kappa0=0.0, kappa1=0.025, kappa2=2.0, # damped investment for stable orbit pi1=0.0, pi2=0.00236, # Nordhaus (2008) damage coefficients tau_T=50.0, lambda_T=1.13, F2x=3.7, CO2_pre=280.0, sigma=0.02, delta_CO2=0.01, n_tax=0.10, ) def damage(T, p): return 1.0 - 1.0 / (1.0 + p["pi1"] * T + p["pi2"] * T**2) def gemmes_ode(state, t, p): omega, lam, d, T, CO2 = state D = damage(T, p) pi = (1.0 - omega - p["r"] * d) * (1.0 - D) phi = p["phi0"] + p["phi1"] / jnp.maximum(1.0 - lam, 0.01)**2 kappa = (p["kappa0"] + p["kappa1"] * jnp.exp(p["kappa2"] * pi)) / p["nu"] F = p["F2x"] * jnp.log(jnp.maximum(CO2, 1.0) / p["CO2_pre"]) / jnp.log(2.0) return jnp.array([ omega * (phi - p["alpha"]), lam * (kappa - p["alpha"] - p["beta"]), kappa - pi - (p["alpha"] + p["beta"]) * d, (F - p["lambda_T"] * T) / p["tau_T"], p["sigma"] * (1.0 - p["n_tax"]) * 100.0 - p["delta_CO2"] * CO2, ]) def simulate(omega0=0.78, lam0=0.95, d0=0.05, T0=0.85, CO2_0=395.0, T_end=80.0, dt=0.02, params=None): p = params or PARAMS state = jnp.array([omega0, lam0, d0, T0, CO2_0]) n = int(round(T_end / dt)) times, traj = [0.0], [state] t = 0.0 for _ in range(n): state = state + gemmes_ode(state, t, p) * dt t += dt times.append(t); traj.append(state) return jnp.array(times), jnp.stack(traj) times, traj = simulate() omega, lam, d, T_anom, CO2 = [traj[:, i] for i in range(5)] print(f"T = {float(times[-1]):.0f} years | {len(times)} steps") print(f"Peak ΔT: {float(T_anom.max()):.2f}°C at t={float(times[int(T_anom.argmax())]):.0f}") print(f"Peak debt: {float(d.max()):.3f} at t={float(times[int(d.argmax())]):.0f}") print(f"Damage at peak ΔT: {float(damage(T_anom.max(), PARAMS))*100:.2f}% of GDP")

T = 80 years  |  4001 steps
Peak ΔT:   1.04°C  at t=24
Peak debt: 0.778     at t=53
Damage at peak ΔT: 0.25% of GDP

In [3]:

fig = plt.figure(figsize=(14, 8))
gs  = gridspec.GridSpec(2, 3, figure=fig, hspace=0.45, wspace=0.35)

ax0 = fig.add_subplot(gs[0, 0])
ax1 = fig.add_subplot(gs[0, 1])
ax2 = fig.add_subplot(gs[0, 2])
ax3 = fig.add_subplot(gs[1, 0])
ax4 = fig.add_subplot(gs[1, 1])
ax5 = fig.add_subplot(gs[1, 2])

t = times

ax0.plot(t, omega, color="steelblue"); ax0.set_title("Wage share ω"); ax0.set_xlabel("Years")
ax1.plot(t, lam,   color="seagreen");  ax1.set_title("Employment λ"); ax1.set_xlabel("Years")
ax2.plot(t, d,     color="crimson")
ax2.axvline(float(times[int(d.argmax())]), ls="--", color="grey", lw=0.8)
ax2.set_title("Debt ratio d"); ax2.set_xlabel("Years")

ax3.plot(t, T_anom, color="darkorange")
ax3.axhline(1.5, ls="--", color="red", lw=0.8, label="Paris 1.5°C")
ax3.axhline(2.0, ls="--", color="darkred", lw=0.8, label="Paris 2.0°C")
ax3.set_title("Temperature anomaly ΔT (°C)"); ax3.set_xlabel("Years"); ax3.legend(fontsize=8)

ax4.plot(t, CO2, color="saddlebrown")
ax4.axhline(280, ls="--", color="grey", lw=0.8, label="pre-industrial")
ax4.set_title("CO₂ (GtC)"); ax4.set_xlabel("Years"); ax4.legend(fontsize=8)

D_path = jnp.array([float(damage(T_anom[i], PARAMS)) * 100 for i in range(len(t))])
ax5.plot(t, D_path, color="purple")
ax5.set_title("Climate damage D(ΔT) (% of GDP)"); ax5.set_xlabel("Years")

fig.suptitle("GEMMES simulation: Keen debt dynamics + climate feedback", fontsize=13)
plt.show()

fig = plt.figure(figsize=(14, 8)) gs = gridspec.GridSpec(2, 3, figure=fig, hspace=0.45, wspace=0.35) ax0 = fig.add_subplot(gs[0, 0]) ax1 = fig.add_subplot(gs[0, 1]) ax2 = fig.add_subplot(gs[0, 2]) ax3 = fig.add_subplot(gs[1, 0]) ax4 = fig.add_subplot(gs[1, 1]) ax5 = fig.add_subplot(gs[1, 2]) t = times ax0.plot(t, omega, color="steelblue"); ax0.set_title("Wage share ω"); ax0.set_xlabel("Years") ax1.plot(t, lam, color="seagreen"); ax1.set_title("Employment λ"); ax1.set_xlabel("Years") ax2.plot(t, d, color="crimson") ax2.axvline(float(times[int(d.argmax())]), ls="--", color="grey", lw=0.8) ax2.set_title("Debt ratio d"); ax2.set_xlabel("Years") ax3.plot(t, T_anom, color="darkorange") ax3.axhline(1.5, ls="--", color="red", lw=0.8, label="Paris 1.5°C") ax3.axhline(2.0, ls="--", color="darkred", lw=0.8, label="Paris 2.0°C") ax3.set_title("Temperature anomaly ΔT (°C)"); ax3.set_xlabel("Years"); ax3.legend(fontsize=8) ax4.plot(t, CO2, color="saddlebrown") ax4.axhline(280, ls="--", color="grey", lw=0.8, label="pre-industrial") ax4.set_title("CO₂ (GtC)"); ax4.set_xlabel("Years"); ax4.legend(fontsize=8) D_path = jnp.array([float(damage(T_anom[i], PARAMS)) * 100 for i in range(len(t))]) ax5.plot(t, D_path, color="purple") ax5.set_title("Climate damage D(ΔT) (% of GDP)"); ax5.set_xlabel("Years") fig.suptitle("GEMMES simulation: Keen debt dynamics + climate feedback", fontsize=13) plt.show()

2. Stranded assets as curvature¶

At the Minsky-climate moment, brown capital is systematically overvalued. The market price of fossil-fuel assets embeds an implicit assumption that they will continue to generate profits — but climate policy and physical damage will strand some of them before they are fully utilised.

The gap between the market price and the climate-adjusted price is a genuine arbitrage: it persists not because markets are irrational but because the risk (policy reversal, physical damage, transition speed) has not yet been priced in.

In gauge-theory language: this is non-zero curvature on the Pacioli manifold. The field strength tensor $F$ measures the arbitrage profit on a round trip through brown and green capital markets.

$$\partial^2(bs) = F \neq 0$$

CurvedBalanceSheet carries this field strength explicitly, making it observable, computable, and part of the calibration target.

In [4]:

# Track the stranded asset premium (curvature ||F||) along the trajectory.
# Premium = climate damage × debt ratio — how much brown capital is overvalued
# relative to its climate-adjusted net present value.

peak_idx = int(d.argmax())
stranded_premiums = []
holonomies        = []

for i in range(0, len(times), 20):
    T_i   = float(T_anom[i])
    d_i   = float(d[i])
    om_i  = float(omega[i])
    D_i   = float(damage(jnp.array(T_i), PARAMS))
    premium = D_i * max(d_i, 0.0) * 50.0   # NPV of stranded loss on brown debt

    positions = jnp.array([
        [ om_i,  0.0,    0.0  ],
        [-om_i, -d_i,    0.2  ],
        [ 0.0,   d_i,   -0.2  ],
        [ 0.0,   premium, 0.0 ],
    ])
    F = jnp.array([0.0, premium, 0.0])

    cbs = CurvedBalanceSheet(
        positions=positions,
        sectors=["workers", "firms", "banks", "climate"],
        instruments=["wages", "brown_capital", "green_capital"],
        curvature=F,
    )
    stranded_premiums.append(float(jnp.linalg.norm(F)))

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].plot(times[::20], stranded_premiums, color="saddlebrown", lw=2)
axes[0].set_xlabel("Years"); axes[0].set_ylabel("||F|| (stranded asset curvature)")
axes[0].set_title("Stranded asset risk (field strength ||F||)\nalong GEMMES trajectory")

axes[1].scatter(d[::20], T_anom[::20], c=times[::20], cmap="plasma", s=8)
axes[1].set_xlabel("Debt ratio d"); axes[1].set_ylabel("Temperature anomaly ΔT (°C)")
axes[1].set_title("Climate-debt phase portrait\ncolour = time")
cb = plt.colorbar(axes[1].collections[0], ax=axes[1])
cb.set_label("Years")

plt.tight_layout()
plt.show()

print(f"Peak stranded asset curvature ||F|| = {max(stranded_premiums):.4f}")
print(f"  (at t ≈ {float(times[::20][stranded_premiums.index(max(stranded_premiums))]):.0f})")

# Track the stranded asset premium (curvature ||F||) along the trajectory. # Premium = climate damage × debt ratio — how much brown capital is overvalued # relative to its climate-adjusted net present value. peak_idx = int(d.argmax()) stranded_premiums = [] holonomies = [] for i in range(0, len(times), 20): T_i = float(T_anom[i]) d_i = float(d[i]) om_i = float(omega[i]) D_i = float(damage(jnp.array(T_i), PARAMS)) premium = D_i * max(d_i, 0.0) * 50.0 # NPV of stranded loss on brown debt positions = jnp.array([ [ om_i, 0.0, 0.0 ], [-om_i, -d_i, 0.2 ], [ 0.0, d_i, -0.2 ], [ 0.0, premium, 0.0 ], ]) F = jnp.array([0.0, premium, 0.0]) cbs = CurvedBalanceSheet( positions=positions, sectors=["workers", "firms", "banks", "climate"], instruments=["wages", "brown_capital", "green_capital"], curvature=F, ) stranded_premiums.append(float(jnp.linalg.norm(F))) fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(times[::20], stranded_premiums, color="saddlebrown", lw=2) axes[0].set_xlabel("Years"); axes[0].set_ylabel("||F|| (stranded asset curvature)") axes[0].set_title("Stranded asset risk (field strength ||F||)\nalong GEMMES trajectory") axes[1].scatter(d[::20], T_anom[::20], c=times[::20], cmap="plasma", s=8) axes[1].set_xlabel("Debt ratio d"); axes[1].set_ylabel("Temperature anomaly ΔT (°C)") axes[1].set_title("Climate-debt phase portrait\ncolour = time") cb = plt.colorbar(axes[1].collections[0], ax=axes[1]) cb.set_label("Years") plt.tight_layout() plt.show() print(f"Peak stranded asset curvature ||F|| = {max(stranded_premiums):.4f}") print(f" (at t ≈ {float(times[::20][stranded_premiums.index(max(stranded_premiums))]):.0f})")

Peak stranded asset curvature ||F|| = 0.0703
  (at t ≈ 52)

3. Green transition as a PCL fold¶

The Keen Minsky moment had two options: invest or deleverage. GEMMES adds a third: invest green — redirect capital from fossil fuel to renewable capacity.

In PCL, this is fold(β, [invest_brown, invest_green, deleverage]):

• β = 0: uniform weights across all three — maximum hedging, no commitment
• β → ∞: winner-take-all — full commitment to the highest-value strategy

The carbon tax shifts the value ranking: without it, brown investment wins on short-run profitability; with it, green investment wins on NPV net of the tax. The β parameter captures how decisively the economy responds to that price signal.

In [5]:

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

wages         = flow("firms", "households", "deposits",  30.0)
invest_brown  = flow("banks", "firms",      "brown_loans", 15.0)
invest_green  = flow("banks", "firms",      "green_loans", 15.0)
deleverage    = flow("firms", "banks",      "brown_loans", 15.0)

# β-sweep: track green_loans balance as β increases
betas = [0.0, 0.5, 1.0, 2.0, 3.0, 5.0, 8.0, 12.0]
brown_out, green_out = [], []

for beta in betas:
    circuit = sequence(wages, fold(beta, [invest_brown, invest_green, deleverage]))
    result  = circuit(bs_pcl)
    assert result.is_consistent(atol=1e-4)
    brown_out.append(float(result.positions[1, 1]))
    green_out.append(float(result.positions[1, 2]))

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(betas, brown_out, "o-", color="saddlebrown", lw=2, label="brown loans (firms)")
ax.plot(betas, green_out, "o-", color="seagreen",    lw=2, label="green loans (firms)")
ax.axhline(-80, ls="--", color="grey", lw=0.8, label="initial brown")
ax.axhline(-20, ls="--", color="lightgreen", lw=0.8, label="initial green")
ax.set_xlabel("β (commitment)"); ax.set_ylabel("Firm loan balance")
ax.set_title("PCL fold(β, [brown, green, deleverage])\nGreen loans rise as commitment increases")
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

bs_pcl = BalanceSheet( positions=jnp.array([ [ 100.0, 0.0, 20.0], [ 0.0, -80.0, -20.0], [-100.0, 80.0, 0.0], ]), sectors=["households", "firms", "banks"], instruments=["deposits", "brown_loans", "green_loans"], ) wages = flow("firms", "households", "deposits", 30.0) invest_brown = flow("banks", "firms", "brown_loans", 15.0) invest_green = flow("banks", "firms", "green_loans", 15.0) deleverage = flow("firms", "banks", "brown_loans", 15.0) # β-sweep: track green_loans balance as β increases betas = [0.0, 0.5, 1.0, 2.0, 3.0, 5.0, 8.0, 12.0] brown_out, green_out = [], [] for beta in betas: circuit = sequence(wages, fold(beta, [invest_brown, invest_green, deleverage])) result = circuit(bs_pcl) assert result.is_consistent(atol=1e-4) brown_out.append(float(result.positions[1, 1])) green_out.append(float(result.positions[1, 2])) fig, ax = plt.subplots(figsize=(8, 4)) ax.plot(betas, brown_out, "o-", color="saddlebrown", lw=2, label="brown loans (firms)") ax.plot(betas, green_out, "o-", color="seagreen", lw=2, label="green loans (firms)") ax.axhline(-80, ls="--", color="grey", lw=0.8, label="initial brown") ax.axhline(-20, ls="--", color="lightgreen", lw=0.8, label="initial green") ax.set_xlabel("β (commitment)"); ax.set_ylabel("Firm loan balance") ax.set_title("PCL fold(β, [brown, green, deleverage])\nGreen loans rise as commitment increases") ax.legend(fontsize=9) plt.tight_layout() plt.show()

4. Carbon tax and TIR routing¶

A carbon tax changes the relative profitability of investment strategies. Without a tax, brown investment wins on short-run profitability. With a tax, the green strategy's NPV improves and brown's worsens.

TIR routing translates these profitability scores into routing weights via the Gibbs distribution. The plot below shows how the routing shifts as we vary the carbon tax rate, at fixed β = 3.

In [6]:

labels = ["brown", "green", "deleverage"]
D_end  = float(damage(T_anom[-1], PARAMS))

# Base scores at end of simulation (climate damage already eroding brown NPV)
base_brown = 0.20 - D_end          # brown eroded by climate damage
base_green = 0.18                  # green unaffected
base_delev = 0.15                  # deleverage always available

# Vary carbon tax rate from 0% to 15%
tax_rates   = [0.00, 0.02, 0.05, 0.08, 0.10, 0.12, 0.15]
w_brown, w_green, w_delev = [], [], []

for tax in tax_rates:
    scores = [base_brown - tax, base_green + tax * 0.6, base_delev]
    w = route(tir_from_scores(labels, scores, beta=3.0))
    w_brown.append(float(w[0]))
    w_green.append(float(w[1]))
    w_delev.append(float(w[2]))

fig, axes = plt.subplots(1, 2, figsize=(13, 4))

# Left: routing weights vs carbon tax rate
axes[0].plot(tax_rates, w_brown, "o-", color="saddlebrown", lw=2, label="brown")
axes[0].plot(tax_rates, w_green, "o-", color="seagreen",    lw=2, label="green")
axes[0].plot(tax_rates, w_delev, "o-", color="steelblue",   lw=2, label="deleverage")
axes[0].set_xlabel("Carbon tax rate"); axes[0].set_ylabel("TIR routing weight (β=3)")
axes[0].set_title("Carbon tax shifts routing from brown to green")
axes[0].legend(); axes[0].set_ylim(0, 1)

# Right: compare no-tax vs 10% tax across β values
beta_sweep = [0.5, 1.0, 2.0, 3.0, 5.0, 8.0]
green_notax, green_tax = [], []
for b in beta_sweep:
    w0 = route(tir_from_scores(labels, [base_brown, base_green, base_delev], beta=b))
    w1 = route(tir_from_scores(labels, [base_brown - 0.10, base_green + 0.06, base_delev], beta=b))
    green_notax.append(float(w0[1]))
    green_tax.append(float(w1[1]))

axes[1].plot(beta_sweep, green_notax, "o--", color="seagreen", lw=2, alpha=0.5, label="no tax")
axes[1].plot(beta_sweep, green_tax,   "o-",  color="seagreen", lw=2,            label="10% carbon tax")
axes[1].fill_between(beta_sweep, green_notax, green_tax, alpha=0.15, color="seagreen")
axes[1].set_xlabel("β"); axes[1].set_ylabel("Green routing weight")
axes[1].set_title("Green routing weight: tax effect amplifies with β\n(more decisive economy responds more to price signal)")
axes[1].legend(); axes[1].set_ylim(0, 1)

plt.tight_layout()
plt.show()

labels = ["brown", "green", "deleverage"] D_end = float(damage(T_anom[-1], PARAMS)) # Base scores at end of simulation (climate damage already eroding brown NPV) base_brown = 0.20 - D_end # brown eroded by climate damage base_green = 0.18 # green unaffected base_delev = 0.15 # deleverage always available # Vary carbon tax rate from 0% to 15% tax_rates = [0.00, 0.02, 0.05, 0.08, 0.10, 0.12, 0.15] w_brown, w_green, w_delev = [], [], [] for tax in tax_rates: scores = [base_brown - tax, base_green + tax * 0.6, base_delev] w = route(tir_from_scores(labels, scores, beta=3.0)) w_brown.append(float(w[0])) w_green.append(float(w[1])) w_delev.append(float(w[2])) fig, axes = plt.subplots(1, 2, figsize=(13, 4)) # Left: routing weights vs carbon tax rate axes[0].plot(tax_rates, w_brown, "o-", color="saddlebrown", lw=2, label="brown") axes[0].plot(tax_rates, w_green, "o-", color="seagreen", lw=2, label="green") axes[0].plot(tax_rates, w_delev, "o-", color="steelblue", lw=2, label="deleverage") axes[0].set_xlabel("Carbon tax rate"); axes[0].set_ylabel("TIR routing weight (β=3)") axes[0].set_title("Carbon tax shifts routing from brown to green") axes[0].legend(); axes[0].set_ylim(0, 1) # Right: compare no-tax vs 10% tax across β values beta_sweep = [0.5, 1.0, 2.0, 3.0, 5.0, 8.0] green_notax, green_tax = [], [] for b in beta_sweep: w0 = route(tir_from_scores(labels, [base_brown, base_green, base_delev], beta=b)) w1 = route(tir_from_scores(labels, [base_brown - 0.10, base_green + 0.06, base_delev], beta=b)) green_notax.append(float(w0[1])) green_tax.append(float(w1[1])) axes[1].plot(beta_sweep, green_notax, "o--", color="seagreen", lw=2, alpha=0.5, label="no tax") axes[1].plot(beta_sweep, green_tax, "o-", color="seagreen", lw=2, label="10% carbon tax") axes[1].fill_between(beta_sweep, green_notax, green_tax, alpha=0.15, color="seagreen") axes[1].set_xlabel("β"); axes[1].set_ylabel("Green routing weight") axes[1].set_title("Green routing weight: tax effect amplifies with β\n(more decisive economy responds more to price signal)") axes[1].legend(); axes[1].set_ylim(0, 1) plt.tight_layout() plt.show()

5. Thermal Shapley: four-player attribution¶

The GEMMES crisis has four potential actors:

Player	Action	Lever
Workers	Moderate wage demands	Reduces debt pressure
Firms	Early green transition	Avoids stranded asset loss
Banks	Tighten brown credit	Limits debt overhang
Policy	Carbon tax	Directly reduces climate damage

Thermal Shapley values attribute the crisis to each player — and identify who has the largest marginal impact (the bottleneck player).

The β parameter in Shapley values has the same interpretation as elsewhere:

• β = 0: classical Shapley (all coalitions equally weighted)
• β → ∞: all weight on the single best coalition — identifies the one player whose action alone would most improve the outcome

In [7]:

peak_damage = float(damage(T_anom.max(), PARAMS))
peak_debt   = float(d.max())

def climate_crisis_value(coalition: frozenset) -> float:
    base = -(peak_damage + 0.1 * peak_debt)
    if 0 in coalition: base += 0.05   # workers
    if 1 in coalition: base += 0.25   # firms
    if 2 in coalition: base += 0.15   # banks
    if 3 in coalition: base += 0.30   # policy
    return base

player_names = ["workers", "firms", "banks", "policy"]
colors_4     = ["steelblue", "crimson", "seagreen", "darkorange"]

beta_sweep_s = [0.0, 0.5, 1.0, 2.0, 4.0, 6.0]
all_phi = []
for b in beta_sweep_s:
    r, _ = full_shapley_analysis(climate_crisis_value, n_players=4, beta=b)
    all_phi.append([float(v) for v in r.phi])
all_phi_T = list(zip(*all_phi))

fig, axes = plt.subplots(1, 2, figsize=(13, 4))

for name, phis, c in zip(player_names, all_phi_T, colors_4):
    axes[0].plot(beta_sweep_s, phis, "o-", color=c, label=name, lw=2)
axes[0].set_xlabel("β"); axes[0].set_ylabel("Shapley value φ")
axes[0].set_title("Thermal Shapley: 4-player climate-finance crisis\nPolicy and firms have largest leverage")
axes[0].legend(fontsize=9)

result_b2, _ = full_shapley_analysis(climate_crisis_value, n_players=4, beta=2.0)
phi_vals = [float(v) for v in result_b2.phi]
bars = axes[1].bar(player_names, phi_vals, color=colors_4)
axes[1].set_ylabel("φ (β=2)")
axes[1].set_title(f"Crisis attribution at β=2\nBottleneck: {player_names[result_b2.bottleneck_player]}")
# Highlight the bottleneck
bars[result_b2.bottleneck_player].set_edgecolor("black")
bars[result_b2.bottleneck_player].set_linewidth(2.5)

plt.tight_layout()
plt.show()

print("Shapley values (β=2):")
for n, v in zip(player_names, phi_vals):
    print(f"  {n:8s}: φ = {v:+.4f}")
print(f"Bottleneck: {player_names[result_b2.bottleneck_player]}")

peak_damage = float(damage(T_anom.max(), PARAMS)) peak_debt = float(d.max()) def climate_crisis_value(coalition: frozenset) -> float: base = -(peak_damage + 0.1 * peak_debt) if 0 in coalition: base += 0.05 # workers if 1 in coalition: base += 0.25 # firms if 2 in coalition: base += 0.15 # banks if 3 in coalition: base += 0.30 # policy return base player_names = ["workers", "firms", "banks", "policy"] colors_4 = ["steelblue", "crimson", "seagreen", "darkorange"] beta_sweep_s = [0.0, 0.5, 1.0, 2.0, 4.0, 6.0] all_phi = [] for b in beta_sweep_s: r, _ = full_shapley_analysis(climate_crisis_value, n_players=4, beta=b) all_phi.append([float(v) for v in r.phi]) all_phi_T = list(zip(*all_phi)) fig, axes = plt.subplots(1, 2, figsize=(13, 4)) for name, phis, c in zip(player_names, all_phi_T, colors_4): axes[0].plot(beta_sweep_s, phis, "o-", color=c, label=name, lw=2) axes[0].set_xlabel("β"); axes[0].set_ylabel("Shapley value φ") axes[0].set_title("Thermal Shapley: 4-player climate-finance crisis\nPolicy and firms have largest leverage") axes[0].legend(fontsize=9) result_b2, _ = full_shapley_analysis(climate_crisis_value, n_players=4, beta=2.0) phi_vals = [float(v) for v in result_b2.phi] bars = axes[1].bar(player_names, phi_vals, color=colors_4) axes[1].set_ylabel("φ (β=2)") axes[1].set_title(f"Crisis attribution at β=2\nBottleneck: {player_names[result_b2.bottleneck_player]}") # Highlight the bottleneck bars[result_b2.bottleneck_player].set_edgecolor("black") bars[result_b2.bottleneck_player].set_linewidth(2.5) plt.tight_layout() plt.show() print("Shapley values (β=2):") for n, v in zip(player_names, phi_vals): print(f" {n:8s}: φ = {v:+.4f}") print(f"Bottleneck: {player_names[result_b2.bottleneck_player]}")

Shapley values (β=2):
  workers : φ = +0.0500
  firms   : φ = +0.2500
  banks   : φ = +0.1500
  policy  : φ = +0.3000
Bottleneck: policy

Summary¶

| Step | Module | What it showed | | --- | --- | --- | | GEMMES ODE | JAX (inline) | 5-variable system; 6-panel trajectory; climate-debt phase portrait | | Stranded assets | CurvedBalanceSheet | Field strength ||F|| peaks at Minsky-climate moment | | Green transition | PCL fold | Three-way β-sweep; green loans rise as commitment increases | | Carbon tax routing | TIR | Tax shifts routing from brown to green; effect amplifies with β | | 4-player Shapley | thermal_shapley | Policy is the bottleneck — largest single lever for mitigation |

Key finding: the carbon tax effect on TIR routing weights amplifies with β. A more decisive economy (higher β) responds more strongly to the price signal — suggesting that institutional commitment mechanisms (green taxonomy, mandatory disclosure) increase the transmission of carbon pricing into investment flows.

---

Previous tutorial: Keen predator-prey

---

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