tegwick
dc4e3adc99
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Assistant: claude-code Assistant-Model: opus Assistant-Process: 14913@bnt-lap001 Assistant-Session: add7ef69-b901-4756-9a67-34cf1125d96f
239 lines
11 KiB
Python
239 lines
11 KiB
Python
"""Calibrate GADM-0.4 starting states to measured data (ASSEMBLY-WP-0013-T02/T03).
|
||
|
||
Targets: calibration_targets.json (ranges per jurisdiction; operators in
|
||
models/Calibration.md). Per jurisdiction:
|
||
|
||
1. Set stocks with targets directly (power channels, E, Sigma) and fit the
|
||
parameters that place them at rest: channel counters (gadm04.calibrate),
|
||
ecological regeneration `ereg` (dE = 0), base labour share `wL0` (claims
|
||
breadth B on target), attention extraction `aext` (attention at rest on
|
||
target, bisection). Where a measured trend of B exists (US), the
|
||
automation strength `auto` is fitted to it. Archetype levers of the nearest archetype
|
||
(models/AdoptionBand.md § 4.4) are applied except `perm`.
|
||
2. Relax wealth, cohort attention and cohort agency to equilibrium.
|
||
3. Adoption: with an anchor's adoption parameters scaled by k, find the model
|
||
time t* where K_Omega's one-year ratio equals the observed ratio, and fit k
|
||
so that K_Omega(t*) hits the level target. Fixed anchors (k = 1) are also
|
||
aligned by growth alone, and their level miss is reported: that is the
|
||
data's verdict on each anchor.
|
||
4. The calibrated start is the full model state at t* with the fitted
|
||
parameters (results/calibrated_starts.json).
|
||
|
||
Residuals are z-scores: (model - centre) / half-width of the target range;
|
||
|z| <= 1 means inside the range.
|
||
|
||
Usage: python3 -B simulations/gadm02/calibrate_data.py
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import json
|
||
import math
|
||
import sys
|
||
from pathlib import Path
|
||
|
||
sys.dont_write_bytecode = True
|
||
|
||
import gadm04 as g # noqa: E402
|
||
import grounding as gr # noqa: E402
|
||
from jurisdictions import ARCHETYPES # noqa: E402
|
||
|
||
HERE = Path(__file__).resolve().parent
|
||
OUT = HERE / "results"
|
||
TARGETS = json.loads((HERE / "calibration_targets.json").read_text())["targets"]
|
||
NEAREST = {"US": "J3", "EU": "J1", "DE": "J1", "AR": "J2", "CN": "J4"}
|
||
FALLBACK = {"DE": "EU"} # DE uses EU targets where it has none
|
||
CHANNEL_KEYS = ("Pm", "Pi", "Pd", "Pmem", "Plin")
|
||
ADOPT = ("aO", "ell", "oseed")
|
||
STEP = 0.1
|
||
BASE_ANCHOR = "B2" # adoption parameters scaled by k start from B2
|
||
|
||
|
||
def targets(jur: str) -> dict:
|
||
t = dict(TARGETS.get(FALLBACK.get(jur, ""), {}))
|
||
t.update(TARGETS[jur])
|
||
return t
|
||
|
||
|
||
def centre(t: dict) -> float:
|
||
return 0.5 * (t["lo"] + t["hi"])
|
||
|
||
|
||
def z(value: float, t: dict) -> float:
|
||
return (value - centre(t)) / (0.5 * (t["hi"] - t["lo"]))
|
||
|
||
|
||
def relax(s: dict, p: dict, keys=gr.EQ_KEYS, steps: int = 4000, h: float = 0.25) -> dict:
|
||
s = dict(s)
|
||
for _ in range(steps):
|
||
d, _a = g.derivatives(s, p)
|
||
for k in keys:
|
||
s[k] = min(1.0, max(0.0, s[k] + h * d[k]))
|
||
return s
|
||
|
||
|
||
def start_state(jur: str, anchor: str, auto: float | None = None) -> tuple[dict, dict]:
|
||
"""(params, init) at rest with the jurisdiction's state targets, before adoption alignment."""
|
||
t = targets(jur)
|
||
vp, vi = gr.variant(gr.BAND[anchor])
|
||
_title, lev, _s0 = ARCHETYPES[NEAREST[jur]]
|
||
p = {**g.PARAMS, **vp, **{k: v for k, v in lev.items() if k != "perm"}}
|
||
if auto is not None:
|
||
p["auto"] = auto
|
||
s = {**g.INIT, **vi}
|
||
for k in CHANNEL_KEYS:
|
||
if k in t:
|
||
s[k] = centre(t[k])
|
||
if "E" in t:
|
||
s["E"] = centre(t["E"])
|
||
if "Sig" in t:
|
||
s["Sig"] = s["Lg"] = centre(t["Sig"])
|
||
p["sig_ref"] = s["Sig"] # BR-009: sentiment anchored at its measured level
|
||
for _ in range(2):
|
||
if "B" in t: # base labour share so that B(t0) is on target
|
||
c = p["dshare"] + p["cp"]
|
||
k0 = c + max(0.0, 1 - c) * p["od"]
|
||
wl = (centre(t["B"]) - k0) / (1 - k0)
|
||
p["wL0"] = max(0.05, min(0.95, wl / (1 - p["auto"] * s["Ko"])))
|
||
d, a = g.derivatives(s, p)
|
||
if "E" in t: # regeneration so that E is at rest
|
||
drain = (p["e1"] * s["Kp"] * a["dem"] * p["thr"] + p["e2"] * s["Ko"])
|
||
p["ereg"] = drain * s["E"] / (1 - s["E"])
|
||
p.update(g.calibrate(p, s))
|
||
s = relax(s, p)
|
||
if "A" in t:
|
||
p["aext"], s = fit_attention(s, p, centre(t["A"]))
|
||
return p, s
|
||
|
||
|
||
def fit_attention(s: dict, p: dict, target: float) -> tuple[float, dict]:
|
||
"""Bisection on the extraction strength aext so that attention at rest equals `target`."""
|
||
lo, hi = 0.05, 3.0
|
||
for _ in range(14):
|
||
mid = 0.5 * (lo + hi)
|
||
s2 = relax(s, {**p, "aext": mid}, steps=1500)
|
||
_d, a = g.derivatives(s2, {**p, "aext": mid})
|
||
lo, hi = (mid, hi) if a["A"] > target else (lo, mid)
|
||
aext = 0.5 * (lo + hi)
|
||
return aext, relax(s, {**p, "aext": aext})
|
||
|
||
|
||
def align(rows: list[dict], ratio: float) -> int | None:
|
||
"""First row after the growth peak whose one-year K_Omega ratio falls to `ratio`."""
|
||
k = int(round(1 / STEP))
|
||
rs = [rows[i + k]["Ko"] / rows[i]["Ko"] for i in range(len(rows) - k)]
|
||
peak = rs.index(max(rs))
|
||
if rs[peak] < ratio:
|
||
return None
|
||
return next((i for i in range(peak, len(rs)) if rs[i] <= ratio), None)
|
||
|
||
|
||
def adoption_path(p: dict, s: dict, kfac: float, anchor_params: dict | None = None) -> list[dict]:
|
||
q = dict(p)
|
||
if anchor_params is not None:
|
||
q.update({k: anchor_params[k] for k in ADOPT if k in anchor_params})
|
||
q["selfenh"] = anchor_params.get("selfenh", 0.0)
|
||
for k in ADOPT:
|
||
q[k] = q[k] * kfac
|
||
return g.simulate(q, s, T=20, record_every=1), q
|
||
|
||
|
||
def fit_k(jur: str, p: dict, s: dict) -> tuple[float, int, list[dict], dict]:
|
||
t = targets(jur)
|
||
ratio, ko_t = centre(t["Ko_ratio_1y"]), t["Ko"]
|
||
best = None
|
||
for i in range(-16, 17): # k in [0.25, 4], log grid
|
||
kf = 2 ** (i / 8)
|
||
rows, q = adoption_path(p, s, kf)
|
||
j = align(rows, ratio)
|
||
if j is None:
|
||
continue
|
||
err = z(rows[j]["Ko"], ko_t) ** 2
|
||
if best is None or err < best[0]:
|
||
best = (err, kf, j, rows, q)
|
||
if best is None:
|
||
raise RuntimeError(f"{jur}: no adoption speed reproduces the observed growth")
|
||
return best[1], best[2], best[3], best[4]
|
||
|
||
|
||
def residuals(jur: str, row: dict, p: dict, rows: list[dict], j: int) -> dict:
|
||
t = targets(jur)
|
||
k = int(round(1 / STEP))
|
||
model = {"Ko": row["Ko"], "Ko_ratio_1y": rows[j + k]["Ko"] / rows[j]["Ko"] if j + k < len(rows) else math.nan,
|
||
"B": row["B"], "A": row["A"], "E": row["E"], "Sig": row["Sig"],
|
||
**{c: row[c] for c in CHANNEL_KEYS}}
|
||
if "B_trend" in t and j + k < len(rows):
|
||
model["B_trend"] = rows[j + k]["B"] - row["B"]
|
||
return {name: (model[name], z(model[name], tt)) for name, tt in t.items() if name in model}
|
||
|
||
|
||
def main() -> None:
|
||
OUT.mkdir(exist_ok=True)
|
||
lines = [f"# {g.MODEL_VERSION}: calibration to measured data", "",
|
||
"Generated by `simulations/gadm02/calibrate_data.py` (ASSEMBLY-WP-0013-T02/T03).",
|
||
"Targets: `calibration_targets.json`; z = (model − centre) / half-width, |z| ≤ 1 inside the",
|
||
f"range. Adoption scaled from anchor {BASE_ANCHOR} by k; t* aligned by observed growth.", ""]
|
||
starts = {}
|
||
lines += ["## Fitted adoption speed and alignment", "",
|
||
"| Jurisdiction | Archetype levers | k | t* (y) | K_Ω(t*) | fitted wL0 | fitted ereg | fitted aext | auto |",
|
||
"|---|---|---|---|---|---|---|---|---|"]
|
||
verdict = ["", "## Fixed anchors aligned by growth: K_Ω level miss", "",
|
||
"| Jurisdiction | " + " | ".join(gr.BAND) + " |", "|---|" + "---|" * len(gr.BAND)]
|
||
resid_lines = []
|
||
resid_all = {}
|
||
for jur in ("US", "EU", "DE", "AR", "CN"):
|
||
t = targets(jur)
|
||
p, s = start_state(jur, BASE_ANCHOR)
|
||
if "Ko_ratio_1y" not in t:
|
||
lines.append(f"| {jur} | {NEAREST[jur]} | — | — | — | {p['wL0']:.3f} | {p['ereg']:.4f} | {p['aext']:.3f} | {p['auto']:.2f} |")
|
||
starts[jur] = {"params": {k: p[k] for k in p if p[k] != g.PARAMS.get(k)}, "init": s, "t_star": None,
|
||
"note": "no measured adoption growth; adoption left at anchor B2 from its start"}
|
||
continue
|
||
kf, j, rows, q = fit_k(jur, p, s)
|
||
if "B_trend" in t: # automation strength so that B's one-year change is on target
|
||
lo, hi = 0.05, 1.2
|
||
for _ in range(8):
|
||
mid = 0.5 * (lo + hi)
|
||
p, s = start_state(jur, BASE_ANCHOR, auto=mid)
|
||
kf, j, rows, q = fit_k(jur, p, s)
|
||
k1 = j + int(round(1 / STEP))
|
||
lo, hi = (lo, mid) if rows[k1]["B"] - rows[j]["B"] < centre(t["B_trend"]) else (mid, hi)
|
||
p, s = start_state(jur, BASE_ANCHOR, auto=0.5 * (lo + hi))
|
||
kf, j, rows, q = fit_k(jur, p, s)
|
||
row = rows[j]
|
||
lines.append(f"| {jur} | {NEAREST[jur]} | {kf:.2f} | {row['t']:.1f} | {row['Ko']:.3f} | {q['wL0']:.3f}"
|
||
f" | {q['ereg']:.4f} | {q['aext']:.3f} | {q['auto']:.2f} |")
|
||
init = {k: row[k] for k in g.STOCKS}
|
||
starts[jur] = {"params": {k: q[k] for k in q if q[k] != g.PARAMS.get(k)}, "init": init,
|
||
"t_star": row["t"], "k": kf}
|
||
res = residuals(jur, row, q, rows, j)
|
||
resid_all[jur] = res
|
||
resid_lines += [f"### {jur}", "", "| Target | Range | Model at t* | z | basis |", "|---|---|---|---|---|"]
|
||
for name, (val, zz) in res.items():
|
||
tt = t[name]
|
||
flag = "" if abs(zz) <= 1 else " ⚠"
|
||
resid_lines.append(f"| {name} | {tt['lo']}–{tt['hi']} | {val:.3f} | {zz:+.2f}{flag} | {tt['basis']} |")
|
||
resid_lines.append("")
|
||
cells = []
|
||
for a in gr.BAND:
|
||
vp, _vi = gr.variant(gr.BAND[a])
|
||
rows_a, _q = adoption_path(p, s, 1.0, anchor_params=vp)
|
||
ja = align(rows_a, centre(t["Ko_ratio_1y"]))
|
||
cells.append("cannot reach growth" if ja is None else
|
||
f"K_Ω {rows_a[ja]['Ko']:.3f}, z {z(rows_a[ja]['Ko'], t['Ko']):+.1f}")
|
||
verdict.append(f"| {jur} | " + " | ".join(cells) + " |")
|
||
n_in = sum(1 for r in resid_all.values() for _v, zz in r.values() if abs(zz) <= 1)
|
||
n_all = sum(len(r) for r in resid_all.values())
|
||
n_meas = sum(1 for jur, r in resid_all.items() for name in r if targets(jur)[name]["basis"] == "measured")
|
||
n_meas_in = sum(1 for jur, r in resid_all.items() for name, (_v, zz) in r.items()
|
||
if targets(jur)[name]["basis"] == "measured" and abs(zz) <= 1)
|
||
lines += ["", f"**Fit summary:** {n_in} of {n_all} targets inside their range; measured targets "
|
||
f"{n_meas_in} of {n_meas}."]
|
||
lines += verdict + ["", "## Residuals at the calibrated start", ""] + resid_lines
|
||
(OUT / "calibration.md").write_text("\n".join(lines))
|
||
(OUT / "calibrated_starts.json").write_text(json.dumps(starts, indent=1, sort_keys=True))
|
||
print("\n".join(lines))
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|