DCT Laboratory — Volume I, Chapter 14

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi']=110

import numpy as np
SEED = 26114
AKK, AKP, BK = 0.90, 0.20, 4.0
APP, BP = 0.80, 2.0
def A_of(g): return np.array([[AKK, AKP],[g, APP]])
def K_ss(g):
    det = (1-AKK)*(1-APP) - AKP*g
    return ((1-APP)*BK + AKP*BP)/det
def rho(g): return float(max(abs(np.linalg.eigvals(A_of(g)))))
def lip(g):   # spectral norm = sqrt(lambda_max(A^T A))
    return float(np.sqrt(max(np.linalg.eigvals(A_of(g).T @ A_of(g)).real)))

def perturb_norm(k=12, g=0.06, d0=(0.0, 5.0)):
    d = np.array(d0)
    for _ in range(k): d = A_of(g) @ d
    return float(np.linalg.norm(d))

def sens_fd(g, h=0.01):
    return (K_ss(g+h) - K_ss(g))/h

def reference_values():
    return {
        "lip_pre": round(lip(0.06), 4),
        "rho_pre": round(rho(0.06), 4),
        "halve_quarters": round(float(np.log(2)/np.log(1/rho(0.06))), 4),
        "K_ss_006": round(K_ss(0.06), 4), "K_ss_007": round(K_ss(0.07), 4),
        "K_ss_008": round(K_ss(0.08), 4), "K_ss_009": round(K_ss(0.09), 4),
        "sens_fd_low":  round(sens_fd(0.06), 4),
        "sens_fd_high": round(sens_fd(0.08), 4),
        "sens_ratio":   round(sens_fd(0.08)/sens_fd(0.06), 4),
        "pert_norm_12": round(perturb_norm(), 4),
    }
if __name__ == "__main__":
    [print(f"{k:16s} {v}") for k,v in reference_values().items()]

g = 0.06
print(f"Lipschitz (spectral norm) L = {lip(g):.4f}   spectral radius ρ = {rho(g):.4f}")
print(f"error half-life: ln2/ln(1/ρ) = {np.log(2)/np.log(1/rho(g)):.4f} quarters")
z = np.array([40.0, 50.0]); ss = np.array([K_ss(g), 0]);  # track K error only for display
errs = []
x = np.array([40.0, 50.0])
target = None
# iterate far to find the fixed point numerically, then measure error decay
for _ in range(4000): x = A_of(g)@x + np.array([BK, BP])
fp = x.copy()
x = np.array([40.0, 50.0])
for k in range(60):
    errs.append(float(np.linalg.norm(x-fp))); x = A_of(g)@x + np.array([BK, BP])
r = errs[40]/errs[39]
print(f"measured per-step error ratio at k=40: {r:.4f}  (→ ρ)")

gs = np.linspace(0.05, 0.092, 200)
Ks = [K_ss(g) for g in gs]
fig, ax = plt.subplots(figsize=(8.2,4.4))
ax.plot(gs, Ks, c="#C8A24B", lw=2.5)
for g in (0.06,0.07,0.08,0.09):
    ax.scatter([g],[K_ss(g)], c="#0B3D2E", zorder=5)
    ax.annotate(f"{K_ss(g):.0f}", (g, K_ss(g)), textcoords="offset points", xytext=(6,6), fontsize=9)
ax.set(xlabel="coupling g", ylabel="steady-state capital K*", 
       title="K*(g) = ((1−a_pp)b_K + a_Kp b_p) / det(g) — seed 26114")
ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"sensitivity (fd, h=0.01): at g=0.06 → {sens_fd(0.06):.1f}   at g=0.08 → {sens_fd(0.08):.1f}   ratio {sens_fd(0.08)/sens_fd(0.06):.1f}")

d = np.array([0.0, 5.0]); path = [d.copy()]
for _ in range(12): d = A_of(0.06)@d; path.append(d.copy())
path = np.array(path)
fig, ax = plt.subplots(figsize=(7.8,4.0))
ax.plot(path[:,0], "o-", c="#0B3D2E", lw=2, ms=4, label="δK (was 0)")
ax.plot(path[:,1], "s-", c="#C8A24B", lw=2, ms=4, label="δp (was 5)")
ax.set(xlabel="quarter", ylabel="perturbation", title="A +5 error in p migrates into K")
ax.legend(frameon=False); ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"‖δ‖ at k=12: {np.linalg.norm(path[12]):.4f}")

ref = reference_values()
expected = {"lip_pre":0.9922,"rho_pre":0.9704,"halve_quarters":23.0814,
 "K_ss_006":150.0,"K_ss_007":200.0,"K_ss_008":300.0,"K_ss_009":600.0,
 "sens_fd_low":5000.0,"sens_fd_high":30000.0,"sens_ratio":6.0,"pert_norm_12":3.0097}
for k,v in expected.items():
    assert abs(ref[k]-v)<5e-4, f"MISMATCH {k}"
    print(f"PASS  {k:16s} {ref[k]}")
print("\nAll checkpoints agree — seed 26114.")