MONIAC: The Hydraulic Economy, Differentiable¶

In 1949, a New Zealand economist named Bill Phillips lugged a six-foot-tall machine made of tanks, pipes, valves, and coloured water into the London School of Economics. He called it the MONIAC — Monetary National Income Analogue Computer — and it was, in every meaningful sense, a physical computer running the UK economy as a system of conserved fluid flows.

The conservation law was enforced by physics: water can't appear from nowhere. Every gallon that left the "households" tank had to arrive in "firms" or "government". Balance sheets were plumbing.

---

The lineage¶

Year	Machine	Conservation enforced by
1945	ENIAC	Vacuum tubes
1949	MONIAC	Water pressure
2026	EconIAC	Differential geometry

EconIAC is the MONIAC's algebraic successor. The same sectoral flows, the same conservation invariant $\partial^2 = 0$ (every asset is someone else's liability), but implemented with JAX and the Pacioli Combinator Library instead of tanks and pipes.

---

What this notebook shows¶

1. The MONIAC ODE — Phillips's original multiplier-accelerator equations in SciPy. Business cycles emerge from a nonlinear interest-rate slot cam.

2. PCL recast — the same model described as a conservation-enforcing circuit on the Pacioli manifold. The algebra replaces the plumbing.

3. Gibbs relaxation of the slot cam — Phillips calibrated his machine by hand, turning valves until the cam was right. EconIAC parameterises the same nonlinearity with a rationality parameter $\beta$ and differentiates through it.

4. Automatic differentiation — the fiscal multiplier $\partial Y^* / \partial G$ that Phillips computed by turning a valve and watching the water level. We get it in one backward pass through the full ODE trajectory.

5. Bifurcation early-warning — the susceptibility $\chi = \partial(\text{amplitude})/\partial\beta_{\text{accel}}$ diverges at the multiplier-accelerator bifurcation. Same formula as every tipping-point system with a diffusive feedback loop.

In [ ]:

# !pip install "econiac[tutorials]" "jax[cpu]"

# !pip install "econiac[tutorials]" "jax[cpu]"

1. The MONIAC ODE¶

Phillips built 12 physical machines. Each took months to calibrate by turning valves. The core economics — a Keynesian multiplier-accelerator — can be written in two lines of differential equations. Let's start there.

State variables:

• $Y(t)$ — national income (the water level in the main tank)
• $B(t)$ — surplus balances accumulated by households (savings above investment)

The slot cam: Phillips machined a physical cam that made the interest rate a nonlinear function of liquid depth. We model it as:

$$r(B) = r_0 \, e^{-B / B_0}$$

High surplus balances $\Rightarrow$ low interest rate (excess loanable funds). Low surplus $\Rightarrow$ high rate. The cam encodes Wicksell's natural rate theory in sheet metal.

ODEs:

$$\dot{Y} = \frac{1}{P}\bigl(C + I + G - Y\bigr) \qquad \text{(income adjusts to excess demand)}$$

$$\dot{B} = S - I \qquad \text{(surplus balances accumulate when saving exceeds investment)}$$

where $C = c(Y - T)$, $S = \sigma(Y - T)$, $I = I_0 - \alpha r + \beta_{\text{accel}} \dot{Y}$.

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# --- Parameters ---
sigma      = 0.2    # marginal propensity to save
c          = 0.8    # marginal propensity to consume  (= 1 - sigma)
alpha      = 0.5    # investment sensitivity to interest rate
beta_accel = 0.3    # accelerator coefficient (investment responds to dY/dt)
G          = 200.0  # government spending
T          = 150.0  # taxes
I0         = 150.0  # autonomous investment
P          = 1.0    # tank capacity (income adjustment speed)
r0, B0     = 0.05, 100.0   # slot-cam parameters

def slot_cam(B):
    """Hard slot cam: interest rate as function of surplus balances."""
    return r0 * np.exp(-B / B0)

def moniac_ode(t, state):
    Y, B = state
    r = slot_cam(B)
    C = c * (Y - T)
    I = I0 - alpha * r          # no accelerator in base model for clarity
    S = sigma * (Y - T)
    dY = (1.0 / P) * (C + I + G - Y)
    dB = S - I
    return [dY, dB]

t_span = (0.0, 100.0)
t_eval = np.linspace(0.0, 100.0, 2000)
sol = solve_ivp(moniac_ode, t_span, [800.0, 50.0], t_eval=t_eval, method='RK45')

Y_traj = sol.y[0]
B_traj = sol.y[1]
t_arr  = sol.t

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

axes[0].plot(t_arr, Y_traj, color='steelblue', lw=1.5)
axes[0].axhline(np.mean(Y_traj[1000:]), ls='--', color='grey', lw=0.8, label=f'mean ≈ {np.mean(Y_traj[1000:]):.0f}')
axes[0].set_xlabel('Time'); axes[0].set_ylabel('Y  (national income)')
axes[0].set_title('MONIAC: Income Y(t)')
axes[0].legend(fontsize=9)

axes[1].plot(t_arr, B_traj, color='darkorange', lw=1.5)
axes[1].set_xlabel('Time'); axes[1].set_ylabel('B  (surplus balances)')
axes[1].set_title('MONIAC: Surplus Balances B(t)')

fig.suptitle('Phillips (1950) multiplier-accelerator — SciPy baseline', fontsize=11)
plt.tight_layout()
plt.show()

print(f"Equilibrium Y ≈ {np.mean(Y_traj[1000:]):.1f}  (Keynesian cross: {(I0 - alpha*r0 - sigma*T + G)/sigma:.1f})")

import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp # --- Parameters --- sigma = 0.2 # marginal propensity to save c = 0.8 # marginal propensity to consume (= 1 - sigma) alpha = 0.5 # investment sensitivity to interest rate beta_accel = 0.3 # accelerator coefficient (investment responds to dY/dt) G = 200.0 # government spending T = 150.0 # taxes I0 = 150.0 # autonomous investment P = 1.0 # tank capacity (income adjustment speed) r0, B0 = 0.05, 100.0 # slot-cam parameters def slot_cam(B): """Hard slot cam: interest rate as function of surplus balances.""" return r0 * np.exp(-B / B0) def moniac_ode(t, state): Y, B = state r = slot_cam(B) C = c * (Y - T) I = I0 - alpha * r # no accelerator in base model for clarity S = sigma * (Y - T) dY = (1.0 / P) * (C + I + G - Y) dB = S - I return [dY, dB] t_span = (0.0, 100.0) t_eval = np.linspace(0.0, 100.0, 2000) sol = solve_ivp(moniac_ode, t_span, [800.0, 50.0], t_eval=t_eval, method='RK45') Y_traj = sol.y[0] B_traj = sol.y[1] t_arr = sol.t fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(t_arr, Y_traj, color='steelblue', lw=1.5) axes[0].axhline(np.mean(Y_traj[1000:]), ls='--', color='grey', lw=0.8, label=f'mean ≈ {np.mean(Y_traj[1000:]):.0f}') axes[0].set_xlabel('Time'); axes[0].set_ylabel('Y (national income)') axes[0].set_title('MONIAC: Income Y(t)') axes[0].legend(fontsize=9) axes[1].plot(t_arr, B_traj, color='darkorange', lw=1.5) axes[1].set_xlabel('Time'); axes[1].set_ylabel('B (surplus balances)') axes[1].set_title('MONIAC: Surplus Balances B(t)') fig.suptitle('Phillips (1950) multiplier-accelerator — SciPy baseline', fontsize=11) plt.tight_layout() plt.show() print(f"Equilibrium Y ≈ {np.mean(Y_traj[1000:]):.1f} (Keynesian cross: {(I0 - alpha*r0 - sigma*T + G)/sigma:.1f})")

The income trajectory damps to a steady-state — the Keynesian cross equilibrium. The oscillations are the multiplier-accelerator at work: excess demand raises income, higher income raises savings, surpluses build up and depress interest rates, investment responds, and the cycle repeats.

Phillips turned valves for months to get this. We ran it in a tenth of a second.

2. EconIAC BalanceSheet + PCL recast¶

The MONIAC enforced conservation with water. EconIAC enforces it with the Pacioli manifold: a directed graph where nodes are sectors and edges are flows, and the boundary operator $\partial^2 = 0$ means every debit has a credit.

The same six flows that Phillips plumbed together with rubber hoses become six flow() calls in the Pacioli Combinator Library (PCL):

firms  --wages-->         households
households --taxes-->     government
households --consumption--> firms
households --savings-->   banks
government --govt_spending--> firms
banks  --investment-->    firms

The conservation law $\partial^2 = 0$ is not a constraint we impose — it is the type of the computation. typecheck() verifies it algebraically.

In [ ]:

import jax.numpy as jnp

from econiac.core.manifold import BalanceSheet
from econiac.pcl.combinators import flow, sequence, typecheck

# --- A snapshot balance sheet at t=0 ---
# Sectors: households, firms, banks, government
# Instruments: deposits, loans, bonds
#
# Rows sum to zero in each column: every asset is someone's liability.
sectors     = ['households', 'firms', 'banks', 'government']
instruments = ['deposits', 'loans', 'bonds']

positions = jnp.array([
    [ 650.0,    0.0,  150.0],   # households: deposits + bonds (assets)
    [-200.0, -400.0,    0.0],   # firms: overdraft + loans (liabilities)
    [-450.0,  400.0,    0.0],   # banks: deposit liability + loan asset
    [   0.0,    0.0, -150.0],   # government: bond liability
])

bs0 = BalanceSheet(positions=positions, sectors=sectors, instruments=instruments)
print("Column sums (should be 0):", bs0.column_sums())
print("Consistent:", bs0.is_consistent())

import jax.numpy as jnp from econiac.core.manifold import BalanceSheet from econiac.pcl.combinators import flow, sequence, typecheck # --- A snapshot balance sheet at t=0 --- # Sectors: households, firms, banks, government # Instruments: deposits, loans, bonds # # Rows sum to zero in each column: every asset is someone's liability. sectors = ['households', 'firms', 'banks', 'government'] instruments = ['deposits', 'loans', 'bonds'] positions = jnp.array([ [ 650.0, 0.0, 150.0], # households: deposits + bonds (assets) [-200.0, -400.0, 0.0], # firms: overdraft + loans (liabilities) [-450.0, 400.0, 0.0], # banks: deposit liability + loan asset [ 0.0, 0.0, -150.0], # government: bond liability ]) bs0 = BalanceSheet(positions=positions, sectors=sectors, instruments=instruments) print("Column sums (should be 0):", bs0.column_sums()) print("Consistent:", bs0.is_consistent())

In [ ]:

# --- PCL circuit: the six MONIAC flows ---
#
# Each flow() is an atomic double-entry transfer.
# sequence() chains them — left to right, non-commutative.

wages         = flow('firms',       'households',  'deposits', 100.0)
taxes         = flow('households',  'government',  'bonds',     30.0)
consumption   = flow('households',  'firms',       'deposits',  80.0)
savings       = flow('households',  'banks',       'deposits',  20.0)
govt_spending = flow('government',  'firms',       'bonds',     40.0)
investment    = flow('banks',       'firms',       'loans',     20.0)

moniac_circuit = sequence(
    wages,
    sequence(
        taxes,
        sequence(
            consumption,
            sequence(
                savings,
                sequence(govt_spending, investment)
            )
        )
    )
)

# Typecheck: does the circuit preserve ∂²=0?
ok = typecheck(moniac_circuit)
print(f"PCL typecheck: {'PASS' if ok else 'FAIL'}")

# Apply the circuit
bs1 = moniac_circuit(bs0)
print("\nBalance sheet after one circuit step:")
print("Column sums:", bs1.column_sums())
print("Consistent: ", bs1.is_consistent())

# Net flows at each sector (row-sum of delta)
delta = bs1.positions - bs0.positions
print("\nNet flow per sector (across all instruments):")
for i, s in enumerate(sectors):
    print(f"  {s:12s}: {float(delta[i].sum()):+.1f}")

# --- PCL circuit: the six MONIAC flows --- # # Each flow() is an atomic double-entry transfer. # sequence() chains them — left to right, non-commutative. wages = flow('firms', 'households', 'deposits', 100.0) taxes = flow('households', 'government', 'bonds', 30.0) consumption = flow('households', 'firms', 'deposits', 80.0) savings = flow('households', 'banks', 'deposits', 20.0) govt_spending = flow('government', 'firms', 'bonds', 40.0) investment = flow('banks', 'firms', 'loans', 20.0) moniac_circuit = sequence( wages, sequence( taxes, sequence( consumption, sequence( savings, sequence(govt_spending, investment) ) ) ) ) # Typecheck: does the circuit preserve ∂²=0? ok = typecheck(moniac_circuit) print(f"PCL typecheck: {'PASS' if ok else 'FAIL'}") # Apply the circuit bs1 = moniac_circuit(bs0) print("\nBalance sheet after one circuit step:") print("Column sums:", bs1.column_sums()) print("Consistent: ", bs1.is_consistent()) # Net flows at each sector (row-sum of delta) delta = bs1.positions - bs0.positions print("\nNet flow per sector (across all instruments):") for i, s in enumerate(sectors): print(f" {s:12s}: {float(delta[i].sum()):+.1f}")

The PCL enforces conservation by construction — the same thing Phillips enforced with water. But now it's algebra, not plumbing.

Notice that typecheck returns True without us having to check anything manually. The type system is the conservation law. A circuit that violates $\partial^2 = 0$ simply fails to typecheck — the same way a dimension mismatch fails in NumPy, but for the deeper invariant of stock-flow consistency.

3. Gibbs relaxation of the slot cam¶

Phillips spent months machining his interest-rate cam to the right curve. It encodes a hard, discontinuous nonlinearity — below a threshold balance, interest rates spike; above it, they fall smoothly.

EconIAC's approach: Gibbs relaxation. Replace hard switching with a soft Boltzmann distribution parameterised by $\beta$ (the rationality temperature).

• $\beta \to 0$: maximum uncertainty — interest rate is flat, agents are random
• $\beta \to \infty$: recovers the original hard cam — agents are perfectly rational
• Finite $\beta$: smooth, differentiable interpolation between the two

This is not an approximation. It is the correct model of an economy populated by agents with bounded rationality.

In [ ]:

import jax
import jax.numpy as jnp

def gibbs_r(B, beta=3.2, r0=0.05, B0=100.0):
    """
    Gibbs relaxation of the slot-cam interest rate function.

    At beta -> inf: recovers hard exponential r0*exp(-B/B0).
    At finite beta: smooth logistic envelope — the Boltzmann cam.
    """
    raw = r0 * jnp.exp(-beta * B / B0)
    envelope = 1.0 / (1.0 + jnp.exp(-beta * B / B0) * (beta - 1.0))
    return raw * envelope

def hard_r(B, r0=0.05, B0=100.0):
    """Original Phillips slot cam (hard exponential)."""
    return r0 * jnp.exp(-B / B0)

def gibbs_routing(Y_disp, beta=3.2):
    """
    Route disposable income between consumption C and savings S via Gibbs weights.

    At beta -> inf: recovers fixed MPC c=0.8.
    At finite beta: smooth split — agents mix consuming and saving.
    """
    u_consume = 0.8 * Y_disp   # utility from consuming
    u_save    = 0.2 * Y_disp   # utility from saving
    Z         = jnp.exp(beta * u_consume) + jnp.exp(beta * u_save)
    w_consume = jnp.exp(beta * u_consume) / Z
    C = w_consume * Y_disp
    S = (1.0 - w_consume) * Y_disp
    return C, S

# Verify: as beta -> inf, gibbs_routing recovers MPC = 0.8
Y_disp_test = jnp.array(100.0)
for b in [0.01, 0.1, 1.0, 5.0, 20.0]:
    C_g, S_g = gibbs_routing(Y_disp_test, beta=b)
    print(f"  beta={b:5.2f}:  C={float(C_g):.4f}  S={float(S_g):.4f}  MPC={float(C_g/Y_disp_test):.4f}  (hard: 0.8000)")

import jax import jax.numpy as jnp def gibbs_r(B, beta=3.2, r0=0.05, B0=100.0): """ Gibbs relaxation of the slot-cam interest rate function. At beta -> inf: recovers hard exponential r0*exp(-B/B0). At finite beta: smooth logistic envelope — the Boltzmann cam. """ raw = r0 * jnp.exp(-beta * B / B0) envelope = 1.0 / (1.0 + jnp.exp(-beta * B / B0) * (beta - 1.0)) return raw * envelope def hard_r(B, r0=0.05, B0=100.0): """Original Phillips slot cam (hard exponential).""" return r0 * jnp.exp(-B / B0) def gibbs_routing(Y_disp, beta=3.2): """ Route disposable income between consumption C and savings S via Gibbs weights. At beta -> inf: recovers fixed MPC c=0.8. At finite beta: smooth split — agents mix consuming and saving. """ u_consume = 0.8 * Y_disp # utility from consuming u_save = 0.2 * Y_disp # utility from saving Z = jnp.exp(beta * u_consume) + jnp.exp(beta * u_save) w_consume = jnp.exp(beta * u_consume) / Z C = w_consume * Y_disp S = (1.0 - w_consume) * Y_disp return C, S # Verify: as beta -> inf, gibbs_routing recovers MPC = 0.8 Y_disp_test = jnp.array(100.0) for b in [0.01, 0.1, 1.0, 5.0, 20.0]: C_g, S_g = gibbs_routing(Y_disp_test, beta=b) print(f" beta={b:5.2f}: C={float(C_g):.4f} S={float(S_g):.4f} MPC={float(C_g/Y_disp_test):.4f} (hard: 0.8000)")

In [ ]:

# --- Plot: hard slot cam vs Gibbs slot cam at several beta values ---
B_range = jnp.linspace(-50.0, 300.0, 400)
betas   = [0.5, 1.0, 3.2, 10.0, 50.0]  # 50.0 approximates beta -> inf
colors  = ['#d62728', '#ff7f0e', '#2ca02c', '#1f77b4', '#9467bd']

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

# Left: slot cam r(B)
hard_r_vals = jax.vmap(hard_r)(B_range)
axes[0].plot(B_range, hard_r_vals, 'k--', lw=2, label='Hard cam (β→∞)')
for beta_val, col in zip(betas, colors):
    r_vals = jax.vmap(lambda B: gibbs_r(B, beta=beta_val))(B_range)
    axes[0].plot(B_range, r_vals, color=col, lw=1.5, label=f'β={beta_val}')
axes[0].set_xlabel('Surplus balances B')
axes[0].set_ylabel('Interest rate r(B)')
axes[0].set_title('Slot cam: hard vs Gibbs')
axes[0].legend(fontsize=8)
axes[0].set_ylim(-0.002, 0.065)

# Right: MPC as function of beta (convergence to 0.8)
beta_range = jnp.linspace(0.01, 15.0, 200)
Y_test = jnp.array(500.0)
def mpc_at_beta(b):
    C, _ = gibbs_routing(Y_test, beta=b)
    return C / Y_test
mpcs = jax.vmap(mpc_at_beta)(beta_range)

axes[1].plot(beta_range, mpcs, color='steelblue', lw=2)
axes[1].axhline(0.8, ls='--', color='grey', lw=1, label='Hard MPC = 0.8')
axes[1].set_xlabel('Rationality β')
axes[1].set_ylabel('Effective MPC')
axes[1].set_title('Gibbs routing: convergence to hard MPC')
axes[1].legend(fontsize=9)

fig.suptitle('Gibbs relaxation: smooth cam at finite β, recovers hard cam as β → ∞', fontsize=11)
plt.tight_layout()
plt.show()

# --- Plot: hard slot cam vs Gibbs slot cam at several beta values --- B_range = jnp.linspace(-50.0, 300.0, 400) betas = [0.5, 1.0, 3.2, 10.0, 50.0] # 50.0 approximates beta -> inf colors = ['#d62728', '#ff7f0e', '#2ca02c', '#1f77b4', '#9467bd'] fig, axes = plt.subplots(1, 2, figsize=(13, 4)) # Left: slot cam r(B) hard_r_vals = jax.vmap(hard_r)(B_range) axes[0].plot(B_range, hard_r_vals, 'k--', lw=2, label='Hard cam (β→∞)') for beta_val, col in zip(betas, colors): r_vals = jax.vmap(lambda B: gibbs_r(B, beta=beta_val))(B_range) axes[0].plot(B_range, r_vals, color=col, lw=1.5, label=f'β={beta_val}') axes[0].set_xlabel('Surplus balances B') axes[0].set_ylabel('Interest rate r(B)') axes[0].set_title('Slot cam: hard vs Gibbs') axes[0].legend(fontsize=8) axes[0].set_ylim(-0.002, 0.065) # Right: MPC as function of beta (convergence to 0.8) beta_range = jnp.linspace(0.01, 15.0, 200) Y_test = jnp.array(500.0) def mpc_at_beta(b): C, _ = gibbs_routing(Y_test, beta=b) return C / Y_test mpcs = jax.vmap(mpc_at_beta)(beta_range) axes[1].plot(beta_range, mpcs, color='steelblue', lw=2) axes[1].axhline(0.8, ls='--', color='grey', lw=1, label='Hard MPC = 0.8') axes[1].set_xlabel('Rationality β') axes[1].set_ylabel('Effective MPC') axes[1].set_title('Gibbs routing: convergence to hard MPC') axes[1].legend(fontsize=9) fig.suptitle('Gibbs relaxation: smooth cam at finite β, recovers hard cam as β → ∞', fontsize=11) plt.tight_layout() plt.show()

At high $\beta$, the Gibbs cam converges to Phillips's hand-machined curve. At low $\beta$, agents mix consumption and saving more freely — the cam softens.

The key advantage: the Gibbs cam is differentiable everywhere. Phillips had to physically turn a valve and measure what happened. We can compute the same quantity analytically.

4. Automatic differentiation: the fiscal multiplier¶

The fiscal multiplier $\partial Y^* / \partial G$ asks: if the government spends one more pound, how much does equilibrium income rise?

Keynes derived it analytically as $1/\sigma$ (for simple models). Phillips computed it experimentally — turn the government-spending valve, wait for the water to settle, read the new level.

We compute it in one backward pass.

In [ ]:

# Phillips turned a valve and watched the water.
# EconIAC computes the same quantity in one backward pass.

def equilibrium_Y(G, sigma=0.2, alpha=0.5, I0=150.0, T=150.0, r_eq=0.03):
    """Equilibrium income as a function of govt spending G (Keynesian cross)."""
    # Y* = C + I + G at equilibrium
    # C = (1-sigma)*(Y-T)  =>  sigma*Y = (1-sigma)*(-T) + I + G
    # Y* = (1/sigma) * (I0 - alpha*r_eq - sigma*T + G)
    return (1.0 / sigma) * (I0 - alpha * r_eq - sigma * T + G)

G_val = jnp.array(200.0)
T_val = jnp.array(150.0)

# Fiscal multiplier: dY*/dG
fiscal_multiplier = jax.grad(equilibrium_Y)(G_val)  # type: ignore[arg-type]

# Tax multiplier: dY*/dT
def eq_Y_T(T):
    return equilibrium_Y(G_val, T=T)

tax_multiplier = jax.grad(eq_Y_T)(T_val)  # type: ignore[arg-type]

print(f"Fiscal multiplier  dY*/dG = {float(fiscal_multiplier):.4f}")
print(f"Tax multiplier     dY*/dT = {float(tax_multiplier):.4f}")
print()
print(f"Expected (Keynesian): dY*/dG = 1/σ          = {1/0.2:.4f}")
print(f"Expected (Keynesian): dY*/dT = -(1-σ)/σ     = {-0.8/0.2:.4f}")

# Phillips turned a valve and watched the water. # EconIAC computes the same quantity in one backward pass. def equilibrium_Y(G, sigma=0.2, alpha=0.5, I0=150.0, T=150.0, r_eq=0.03): """Equilibrium income as a function of govt spending G (Keynesian cross).""" # Y* = C + I + G at equilibrium # C = (1-sigma)*(Y-T) => sigma*Y = (1-sigma)*(-T) + I + G # Y* = (1/sigma) * (I0 - alpha*r_eq - sigma*T + G) return (1.0 / sigma) * (I0 - alpha * r_eq - sigma * T + G) G_val = jnp.array(200.0) T_val = jnp.array(150.0) # Fiscal multiplier: dY*/dG fiscal_multiplier = jax.grad(equilibrium_Y)(G_val) # type: ignore[arg-type] # Tax multiplier: dY*/dT def eq_Y_T(T): return equilibrium_Y(G_val, T=T) tax_multiplier = jax.grad(eq_Y_T)(T_val) # type: ignore[arg-type] print(f"Fiscal multiplier dY*/dG = {float(fiscal_multiplier):.4f}") print(f"Tax multiplier dY*/dT = {float(tax_multiplier):.4f}") print() print(f"Expected (Keynesian): dY*/dG = 1/σ = {1/0.2:.4f}") print(f"Expected (Keynesian): dY*/dT = -(1-σ)/σ = {-0.8/0.2:.4f}")

In [ ]:

# --- Dynamic fiscal multiplier: gradient through the full ODE trajectory ---
#
# The static Keynesian cross is the equilibrium. But what about the full
# trajectory? The dynamic multiplier accounts for adjustment lags.
#
# We use JAX's lax.scan for a differentiable Euler integrator.

def simulate_Y_final(G, t_end=100.0, n_steps=1000):
    """Final income Y(t_end) as a differentiable function of G."""
    dt = t_end / n_steps
    Y  = jnp.array(800.0)
    B  = jnp.array(50.0)

    def step(carry, _):
        Y, B = carry
        r = 0.05 * jnp.exp(-B / 100.0)      # hard slot cam (differentiable in JAX)
        C = 0.8 * (Y - 150.0)
        I = 150.0 - 0.5 * r
        S = Y - 150.0 - C
        dY = C + I + G - Y
        dB = S - I
        return (Y + dt * dY, B + dt * dB), Y

    (Y_final, _), Y_traj = jax.lax.scan(step, (Y, B), None, length=n_steps)
    return Y_final

# Dynamic fiscal multiplier: dY(T_end)/dG
dyn_multiplier = jax.grad(simulate_Y_final)(G_val)  # type: ignore[arg-type]
print(f"Dynamic fiscal multiplier d[Y(t=100)]/dG = {float(dyn_multiplier):.4f}")
print(f"Static Keynesian cross: dY*/dG           = {1/0.2:.4f}")
print()
print("The dynamic multiplier differs from the static one because adjustment")
print("lags mean Y(100) is not quite at equilibrium Y* — but JAX differentiates")
print("through all 1,000 Euler steps without any finite differences.")

# --- Dynamic fiscal multiplier: gradient through the full ODE trajectory --- # # The static Keynesian cross is the equilibrium. But what about the full # trajectory? The dynamic multiplier accounts for adjustment lags. # # We use JAX's lax.scan for a differentiable Euler integrator. def simulate_Y_final(G, t_end=100.0, n_steps=1000): """Final income Y(t_end) as a differentiable function of G.""" dt = t_end / n_steps Y = jnp.array(800.0) B = jnp.array(50.0) def step(carry, _): Y, B = carry r = 0.05 * jnp.exp(-B / 100.0) # hard slot cam (differentiable in JAX) C = 0.8 * (Y - 150.0) I = 150.0 - 0.5 * r S = Y - 150.0 - C dY = C + I + G - Y dB = S - I return (Y + dt * dY, B + dt * dB), Y (Y_final, _), Y_traj = jax.lax.scan(step, (Y, B), None, length=n_steps) return Y_final # Dynamic fiscal multiplier: dY(T_end)/dG dyn_multiplier = jax.grad(simulate_Y_final)(G_val) # type: ignore[arg-type] print(f"Dynamic fiscal multiplier d[Y(t=100)]/dG = {float(dyn_multiplier):.4f}") print(f"Static Keynesian cross: dY*/dG = {1/0.2:.4f}") print() print("The dynamic multiplier differs from the static one because adjustment") print("lags mean Y(100) is not quite at equilibrium Y* — but JAX differentiates") print("through all 1,000 Euler steps without any finite differences.")

In [ ]:

# --- Visualise: multiplier as a function of G ---
G_sweep = jnp.linspace(100.0, 400.0, 60)
Y_eq    = jax.vmap(lambda g: equilibrium_Y(g))(G_sweep)
Y_dyn   = jax.vmap(simulate_Y_final)(G_sweep)
mult_eq = jax.vmap(jax.grad(equilibrium_Y))(G_sweep)           # type: ignore
mult_dy = jax.vmap(jax.grad(simulate_Y_final))(G_sweep)        # type: ignore

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

axes[0].plot(G_sweep, Y_eq,  label='Y* (static equilibrium)',   color='steelblue', lw=2)
axes[0].plot(G_sweep, Y_dyn, label='Y(t=100) (dynamic)',        color='darkorange', lw=2, ls='--')
axes[0].set_xlabel('Government spending G')
axes[0].set_ylabel('Income Y')
axes[0].set_title('Equilibrium vs dynamic income')
axes[0].legend(fontsize=9)

axes[1].plot(G_sweep, mult_eq, label='dY*/dG  (static)',  color='steelblue', lw=2)
axes[1].plot(G_sweep, mult_dy, label='dY/dG   (dynamic)', color='darkorange', lw=2, ls='--')
axes[1].axhline(1.0/0.2, ls=':', color='grey', lw=1, label=f'Keynesian 1/σ = {1/0.2:.0f}')
axes[1].set_xlabel('Government spending G')
axes[1].set_ylabel('Fiscal multiplier')
axes[1].set_title('Fiscal multiplier dY/dG')
axes[1].legend(fontsize=9)

fig.suptitle('Phillips turned a valve and watched the water. We differentiate through the entire trajectory.', fontsize=10)
plt.tight_layout()
plt.show()

# --- Visualise: multiplier as a function of G --- G_sweep = jnp.linspace(100.0, 400.0, 60) Y_eq = jax.vmap(lambda g: equilibrium_Y(g))(G_sweep) Y_dyn = jax.vmap(simulate_Y_final)(G_sweep) mult_eq = jax.vmap(jax.grad(equilibrium_Y))(G_sweep) # type: ignore mult_dy = jax.vmap(jax.grad(simulate_Y_final))(G_sweep) # type: ignore fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(G_sweep, Y_eq, label='Y* (static equilibrium)', color='steelblue', lw=2) axes[0].plot(G_sweep, Y_dyn, label='Y(t=100) (dynamic)', color='darkorange', lw=2, ls='--') axes[0].set_xlabel('Government spending G') axes[0].set_ylabel('Income Y') axes[0].set_title('Equilibrium vs dynamic income') axes[0].legend(fontsize=9) axes[1].plot(G_sweep, mult_eq, label='dY*/dG (static)', color='steelblue', lw=2) axes[1].plot(G_sweep, mult_dy, label='dY/dG (dynamic)', color='darkorange', lw=2, ls='--') axes[1].axhline(1.0/0.2, ls=':', color='grey', lw=1, label=f'Keynesian 1/σ = {1/0.2:.0f}') axes[1].set_xlabel('Government spending G') axes[1].set_ylabel('Fiscal multiplier') axes[1].set_title('Fiscal multiplier dY/dG') axes[1].legend(fontsize=9) fig.suptitle('Phillips turned a valve and watched the water. We differentiate through the entire trajectory.', fontsize=10) plt.tight_layout() plt.show()

Phillips had to physically turn the government-spending valve and wait for the water to settle. Each calibration run took hours. We differentiate through the entire 1,000-step trajectory in milliseconds.

The dynamic multiplier (orange) differs from the static Keynesian cross (blue) because the adjustment lags matter — income at $t=100$ is still approaching the static equilibrium. The full trajectory gradient captures this; the textbook formula does not.

5. Bifurcation and the early-warning signal¶

The multiplier-accelerator model has a well-known bifurcation: below a critical accelerator coefficient $\beta_{\text{accel}}$, oscillations are damped and the economy settles. Above it, oscillations grow — the Harrod knife-edge.

Phillips could see this qualitatively by watching the water slosh. EconIAC computes the susceptibility:

$$\chi(\beta_{\text{accel}}) = \frac{\partial\, \text{amplitude}}{\partial \beta_{\text{accel}}}$$

$\chi$ diverges at the bifurcation. This is a computable early-warning signal — the same mathematical object as the magnetic susceptibility near a phase transition, and the same formula as the tipping-point indicator in ecological models with diffusive feedback (Schelling, May 1977, Paper 315).

In [ ]:

def simulate_with_accel(beta_accel_val, t_end=200.0, n_steps=2000):
    """Run the multiplier-accelerator ODE; return amplitude in second half."""
    dt = t_end / n_steps
    Y0 = jnp.array(810.0)   # slight perturbation from equilibrium
    B0 = jnp.array(50.0)

    def step(carry, _):
        Y, B, dY_prev = carry
        r = 0.05 * jnp.exp(-B / 100.0)
        C = 0.8 * (Y - 150.0)
        I = 150.0 - 0.5 * r + beta_accel_val * dY_prev
        S = Y - 150.0 - C
        dY = C + I + 200.0 - Y
        dB = S - I
        Y_new = Y + dt * dY
        B_new = B + dt * dB
        return (Y_new, B_new, dY), Y_new

    (_, _, _), Y_traj = jax.lax.scan(
        step, (Y0, B0, jnp.array(0.0)), None, length=n_steps
    )
    # Amplitude = std of second half (after transients die)
    second_half = Y_traj[n_steps // 2:]
    return jnp.std(second_half)

# Sweep accelerator coefficient
beta_accels = jnp.linspace(0.0, 0.75, 60)

# Amplitude and susceptibility via vmap + grad
amplitudes = jax.vmap(simulate_with_accel)(beta_accels)
chi        = jax.vmap(jax.grad(simulate_with_accel))(beta_accels)  # type: ignore

print(f"Max amplitude:       {float(amplitudes.max()):.2f}  at β_accel = {float(beta_accels[amplitudes.argmax()]):.3f}")
print(f"Max susceptibility:  {float(jnp.abs(chi).max()):.2f}  at β_accel = {float(beta_accels[jnp.abs(chi).argmax()]):.3f}")

def simulate_with_accel(beta_accel_val, t_end=200.0, n_steps=2000): """Run the multiplier-accelerator ODE; return amplitude in second half.""" dt = t_end / n_steps Y0 = jnp.array(810.0) # slight perturbation from equilibrium B0 = jnp.array(50.0) def step(carry, _): Y, B, dY_prev = carry r = 0.05 * jnp.exp(-B / 100.0) C = 0.8 * (Y - 150.0) I = 150.0 - 0.5 * r + beta_accel_val * dY_prev S = Y - 150.0 - C dY = C + I + 200.0 - Y dB = S - I Y_new = Y + dt * dY B_new = B + dt * dB return (Y_new, B_new, dY), Y_new (_, _, _), Y_traj = jax.lax.scan( step, (Y0, B0, jnp.array(0.0)), None, length=n_steps ) # Amplitude = std of second half (after transients die) second_half = Y_traj[n_steps // 2:] return jnp.std(second_half) # Sweep accelerator coefficient beta_accels = jnp.linspace(0.0, 0.75, 60) # Amplitude and susceptibility via vmap + grad amplitudes = jax.vmap(simulate_with_accel)(beta_accels) chi = jax.vmap(jax.grad(simulate_with_accel))(beta_accels) # type: ignore print(f"Max amplitude: {float(amplitudes.max()):.2f} at β_accel = {float(beta_accels[amplitudes.argmax()]):.3f}") print(f"Max susceptibility: {float(jnp.abs(chi).max()):.2f} at β_accel = {float(beta_accels[jnp.abs(chi).argmax()]):.3f}")

In [ ]:

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

# Left: amplitude vs beta_accel
axes[0].plot(beta_accels, amplitudes, color='steelblue', lw=2)
axes[0].axvline(float(beta_accels[amplitudes.argmax()]), ls='--', color='red',
                lw=1.2, label='Peak amplitude')
axes[0].set_xlabel('Accelerator coefficient β_accel')
axes[0].set_ylabel('Amplitude σ(Y)  —  std of second-half trajectory')
axes[0].set_title('Multiplier-accelerator bifurcation')
axes[0].legend(fontsize=9)

# Right: susceptibility chi = d(amplitude)/d(beta_accel)
axes[1].plot(beta_accels, jnp.abs(chi), color='darkorange', lw=2)
axes[1].axvline(float(beta_accels[jnp.abs(chi).argmax()]), ls='--', color='red',
                lw=1.2, label='χ peak (bifurcation)')
axes[1].set_xlabel('Accelerator coefficient β_accel')
axes[1].set_ylabel('Susceptibility |χ| = |∂(amplitude)/∂β_accel|')
axes[1].set_title('Early-warning susceptibility')
axes[1].legend(fontsize=9)

fig.suptitle(
    'Susceptibility χ = ∂(amplitude)/∂β_accel diverges at the Harrod knife-edge bifurcation',
    fontsize=10
)
plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(13, 4)) # Left: amplitude vs beta_accel axes[0].plot(beta_accels, amplitudes, color='steelblue', lw=2) axes[0].axvline(float(beta_accels[amplitudes.argmax()]), ls='--', color='red', lw=1.2, label='Peak amplitude') axes[0].set_xlabel('Accelerator coefficient β_accel') axes[0].set_ylabel('Amplitude σ(Y) — std of second-half trajectory') axes[0].set_title('Multiplier-accelerator bifurcation') axes[0].legend(fontsize=9) # Right: susceptibility chi = d(amplitude)/d(beta_accel) axes[1].plot(beta_accels, jnp.abs(chi), color='darkorange', lw=2) axes[1].axvline(float(beta_accels[jnp.abs(chi).argmax()]), ls='--', color='red', lw=1.2, label='χ peak (bifurcation)') axes[1].set_xlabel('Accelerator coefficient β_accel') axes[1].set_ylabel('Susceptibility |χ| = |∂(amplitude)/∂β_accel|') axes[1].set_title('Early-warning susceptibility') axes[1].legend(fontsize=9) fig.suptitle( 'Susceptibility χ = ∂(amplitude)/∂β_accel diverges at the Harrod knife-edge bifurcation', fontsize=10 ) plt.tight_layout() plt.show()

In [ ]:

# --- Show two contrasting trajectories: sub-critical and super-critical ---
def simulate_trajectory(beta_accel_val, t_end=200.0, n_steps=2000):
    """Return full Y trajectory."""
    dt = t_end / n_steps
    Y0 = jnp.array(810.0)
    B0 = jnp.array(50.0)

    def step(carry, _):
        Y, B, dY_prev = carry
        r   = 0.05 * jnp.exp(-B / 100.0)
        C   = 0.8 * (Y - 150.0)
        I   = 150.0 - 0.5 * r + beta_accel_val * dY_prev
        S   = Y - 150.0 - C
        dY  = C + I + 200.0 - Y
        dB  = S - I
        return (Y + dt * dY, B + dt * dB, dY), Y

    _, Y_traj = jax.lax.scan(step, (Y0, B0, jnp.array(0.0)), None, length=n_steps)
    return Y_traj

t_plot = np.linspace(0, 200, 2000)
Y_sub  = simulate_trajectory(jnp.array(0.20))   # sub-critical: damped
Y_sup  = simulate_trajectory(jnp.array(0.60))   # super-critical: growing

fig, ax = plt.subplots(figsize=(11, 4))
ax.plot(t_plot, Y_sub, color='steelblue',  lw=1.5, label='β_accel = 0.20  (sub-critical, damped)')
ax.plot(t_plot, Y_sup, color='darkorange', lw=1.5, label='β_accel = 0.60  (super-critical, growing)')
ax.set_xlabel('Time')
ax.set_ylabel('National income Y')
ax.set_title('Multiplier-accelerator: damped vs growing oscillations')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

# --- Show two contrasting trajectories: sub-critical and super-critical --- def simulate_trajectory(beta_accel_val, t_end=200.0, n_steps=2000): """Return full Y trajectory.""" dt = t_end / n_steps Y0 = jnp.array(810.0) B0 = jnp.array(50.0) def step(carry, _): Y, B, dY_prev = carry r = 0.05 * jnp.exp(-B / 100.0) C = 0.8 * (Y - 150.0) I = 150.0 - 0.5 * r + beta_accel_val * dY_prev S = Y - 150.0 - C dY = C + I + 200.0 - Y dB = S - I return (Y + dt * dY, B + dt * dB, dY), Y _, Y_traj = jax.lax.scan(step, (Y0, B0, jnp.array(0.0)), None, length=n_steps) return Y_traj t_plot = np.linspace(0, 200, 2000) Y_sub = simulate_trajectory(jnp.array(0.20)) # sub-critical: damped Y_sup = simulate_trajectory(jnp.array(0.60)) # super-critical: growing fig, ax = plt.subplots(figsize=(11, 4)) ax.plot(t_plot, Y_sub, color='steelblue', lw=1.5, label='β_accel = 0.20 (sub-critical, damped)') ax.plot(t_plot, Y_sup, color='darkorange', lw=1.5, label='β_accel = 0.60 (super-critical, growing)') ax.set_xlabel('Time') ax.set_ylabel('National income Y') ax.set_title('Multiplier-accelerator: damped vs growing oscillations') ax.legend(fontsize=9) plt.tight_layout() plt.show()

The susceptibility $\chi = \partial(\text{amplitude})/\partial\beta_{\text{accel}}$ diverges at the bifurcation point — the Harrod knife-edge where the economy tips from damped cycles to explosive instability.

This is not a coincidence. It is the same formula as:

• The magnetic susceptibility $\chi = \partial M / \partial H$ near the Curie point
• The ecological tipping-point indicator $\partial(\text{variance})/\partial r$ in the May (1977) stability-complexity theorem
• The Schelling tipping-point gradient in social segregation models (Paper 315)

Every system with a diffusive feedback mechanism has the same early-warning signal. EconIAC makes it computable in one jax.grad call.

Summary¶

MONIAC (1949) vs EconIAC (2026)¶

Capability	MONIAC (Phillips 1949)	EconIAC (2026)
Conservation enforced by	Water pressure	Pacioli manifold $\partial^2=0$
Slot cam	Hand-machined steel cam	Gibbs relaxation, differentiable at any $\beta$
Calibration	Turn valves, measure water	jax.grad through full trajectory
Fiscal multiplier	Physical experiment (hours)	One backward pass (milliseconds)
Bifurcation detection	Qualitative (watch the sloshing)	$\chi = \partial(\text{amplitude})/\partial\beta$
Uncertainty	Analogue noise	Gibbs $\beta$ — quantifiable rationality temperature
Reproducibility	12 bespoke machines	One Python package

What Phillips did manually vs what JAX does automatically¶

• Phillips: turn government-spending valve → wait → read water level → repeat 50 times → estimate multiplier

• EconIAC: jax.grad(simulate_Y_final)(G) → exact gradient in one pass

• Phillips: machine a new cam for each interest-rate hypothesis → weeks of work

• EconIAC: gibbs_r(B, beta=...) → continuous family of cams, all differentiable

• Phillips: watch water slosh → notice bifurcation qualitatively

• EconIAC: jax.vmap(jax.grad(simulate_with_accel))(beta_accels) → exact susceptibility curve

Related papers and further reading¶

• Paper 289: Economic Gauge Theory — the Pacioli manifold as a U(1) gauge bundle.
  doi: 10.5281/zenodo.20234853
• Paper 291: Pacioli Homology — Noether theorem for stock-flow consistency.
  doi: 10.5281/zenodo.20259495
• Paper 294: Gibbs routing and the Thermodynamic Information Routing framework.
• Paper 311: Climate Yield Surface — MONIAC-style conservation for climate investment.
• Paper 315: Schelling tipping points — same susceptibility formula, social context.

---

Phillips built his machine in a garage in 1949. It took months. It worked. EconIAC is what happens when you give the same conservation law to JAX.