DCT Laboratory — Volume II, Chapter 15

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

import numpy as np
SEED = 26215
A, U, K0, T = 0.9, 1.0, 10.0, 12
SHOCK_T, SHOCK = 8, -6.0
rng = np.random.default_rng([2, SEED])
NOISE = np.round(rng.normal(0, 0.4, T), 4)   # frozen observation noise = canonical
def truth_path():
    K, path = K0, [K0]
    for k in range(T):
        K = A*K + U + (SHOCK if k == SHOCK_T else 0.0)
        path.append(K)
    return np.array(path)
def twin_path(g):
    """twin_{k+1} = (1-g)*(A*twin_k + U) + g*y_{k+1}, y = truth + noise."""
    tr = truth_path()
    tw = [K0]
    for k in range(T):
        y = tr[k+1] + NOISE[k]
        tw.append((1-g)*(A*tw[k] + U) + g*y)
    return np.array(tw)
def rmse(g):
    return float(np.sqrt(np.mean((twin_path(g)-truth_path())**2)))
# decision layer: switch family from estimated terminal state
def J_switch(m, K): return 3.0*(6-m)*np.sqrt(K + 3.0*m)
def m_star(K): return int(np.argmax([J_switch(m, K) for m in range(7)]))
def decision_block():
    Kt = truth_path()[-1]
    Ks, Ko = twin_path(0.5)[-1], twin_path(0.0)[-1]
    ms, mo = m_star(Ks), m_star(Ko)
    return Kt, Ks, Ko, ms, mo, J_switch(ms, Kt) - J_switch(mo, Kt)
# autonomous loop: K_{k+1} = A*K + u, u = 0.1*Ktar + c*(Ktar - K)  =>  gap factor |A - c|
KTAR = 8.0
def gap_after(n, c, K0_=10.0):
    gap = K0_ - KTAR
    return abs(gap*(A-c)**n)
def c_max_stable(step=0.2):
    cs = [round(step*i,1) for i in range(0, 13)]
    return max(c for c in cs if abs(A-c) < 1.0)

def reference_values():
    Kt, Ks, Ko, ms, mo, reg = decision_block()
    return {
        "rmse_open": round(rmse(0.0),4),
        "rmse_sync": round(rmse(0.5),4),
        "sync_advantage": round(rmse(0.0)/rmse(0.5),4),
        "contraction_factor": round((1-0.5)*A,4),
        "K_true_T": round(Kt,4), "K_sync_T": round(Ks,4), "K_open_T": round(Ko,4),
        "mstar_sync": ms, "mstar_open": mo,
        "regret_open_twin": round(reg,4),
        "c_fastest": 0.9, "factor_at_09": 0.0,
        "c_max_stable": c_max_stable(),
        "gap3_c05": round(gap_after(3, 0.5),4),
    }
if __name__ == "__main__":
    print("NOISE:", list(NOISE))
    print("truth:", [round(x,4) for x in truth_path()])
    [print(f"{k:20s} {v}") for k,v in reference_values().items()]

tr = truth_path()
fig, ax = plt.subplots(figsize=(7.8,4.2))
ax.plot(tr, "k--", lw=2, label="enterprise truth (shock −6 at k=8)")
ax.plot(twin_path(0.0), "s-", c="#8A8F8B", lw=1.8, ms=4, label=f"pure-model twin g=0 → RMSE {rmse(0.0):.2f}")
ax.plot(twin_path(0.5), "o-", c="#C8A24B", lw=2.2, ms=4, label=f"synchronized twin g=0.5 → RMSE {rmse(0.5):.2f}")
ax.axvline(SHOCK_T+1, c="#B0532F", ls=":", lw=1.4)
ax.set(xlabel="step", ylabel="capital K", title="The twin that looks vs the twin that remembers — seed 26215")
ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"contraction factor (1-g)*A = {(1-0.5)*A:.4f}   sync advantage: {rmse(0.0)/rmse(0.5):.4f}x")

Kt, Ks, Ko, ms, mo, reg = decision_block()
ms_grid = range(7)
fig, ax = plt.subplots(figsize=(7.8,4.0))
ax.plot(list(ms_grid), [J_switch(m, Kt) for m in ms_grid], "o-", c="#0B3D2E", lw=2.2, label=f"J(m) at TRUE K = {Kt:.2f}")
ax.scatter([ms],[J_switch(ms,Kt)], s=140, c="#C8A24B", zorder=5, label=f"sync twin's pick m={ms}")
ax.scatter([mo],[J_switch(mo,Kt)], s=140, c="#B0532F", marker="X", zorder=5, label=f"drifted twin's pick m={mo}")
ax.set(xlabel="m (build periods before harvest)", ylabel="value at true state",
       title=f"Same optimizer, different twins: regret {reg:.4f} — seed 26215")
ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"K true {Kt:.4f} | sync est {Ks:.4f} -> m*={ms} | open est {Ko:.4f} -> m*={mo} | regret {reg:.4f}")

cs = [round(0.2*i,1) for i in range(13)]
facs = [abs(A-c) for c in cs]
fig, ax = plt.subplots(figsize=(7.8,4.0))
ax.plot(cs, facs, "o-", c="#C8A24B", lw=2.2, ms=5)
ax.axhline(1.0, c="#B0532F", ls=":", lw=1.5, label="stability boundary")
ax.scatter([0.9],[0.0], s=130, c="#0B3D2E", zorder=5, label="deadbeat c = 0.9")
ax.set(xlabel="feedback gain c", ylabel="gap factor |0.9 − c|", title="Stable below the line; oversteering above — seed 26215")
ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"fastest c: 0.9 (factor 0.0)   last stable on grid: {c_max_stable()}   gap after 3 steps at c=0.5: {gap_after(3,0.5):.4f}")

ref = reference_values()
expected = {"rmse_open":2.8811,"rmse_sync":0.9313,"sync_advantage":3.0938,
 "contraction_factor":0.45,"K_true_T":5.626,"K_sync_T":5.8598,"K_open_T":10.0,
 "mstar_sync":1,"mstar_open":0,"regret_open_twin":1.3605,
 "c_fastest":0.9,"factor_at_09":0.0,"c_max_stable":1.8,"gap3_c05":0.128}
for k,v in expected.items():
    assert abs(ref[k]-v)<5e-4, f"MISMATCH {k}"
    print(f"PASS  {k:20s} {ref[k]}")
print("\nAll checkpoints agree — seed 26215.")