1import json
 2import math
 3import numpy as np
 4from scipy.stats import skew, kurtosis
 5
 6SEED = 1499
 7N = 200000
 8EPS = np.array([0.01, 0.02, 0.04, 0.08, 0.16, 0.32])
 9ALPHAS = [0.5, 1.0, 1.5]
10
11def slope(x, y):
12    return float(np.polyfit(np.log(x), np.log(np.maximum(y, 1e-30)), 1)[0])
13
14def sigma2(e, a, c=1.): return c * e**(2-a)/(2-a)
15def rho3(e, a, c=1.): return c * e**(3-a)/(3-a)
16def jump_rate(e, a, c=1.): return c * e**(-a)/a
17
18def main():
19    out = {"N": N, "eps": EPS.tolist(), "results": {}}
20    for a in ALPHAS:
21        rng = np.random.default_rng(SEED + round(100*a))
22        rows = []
23        for e in EPS:
24            s2, r3 = sigma2(e,a), rho3(e,a)
25            # Directly sample the proposed Gaussian replacement. This is a
26            # cheap implementation of the small-jump line in a neural SDE.
27            g = rng.normal(0., math.sqrt(s2), N)
28            z = g / math.sqrt(s2)
29            rows.append({
30              "epsilon": float(e),
31              "sigma2_theory": s2,
32              "gaussian_variance_mc": float(np.var(g)),
33              "variance_relative_error": float(abs(np.var(g)-s2)/s2),
34              # These are the paper's quantitative rate proxies:
35              # naive residual scale sigma; compensated Berry-Esseen/W1 scale
36              # rho/sigma^2; standardized residual normality scale rho/sigma^3.
37              "naive_scale_sigma": math.sqrt(s2),
38              "compensated_w1_proxy_rho_over_sigma2": r3/s2,
39              "normality_proxy_rho_over_sigma3": r3/(s2**1.5),
40              "normal_sample_abs_skew": float(abs(skew(z))),
41              "normal_sample_excess_kurtosis": float(kurtosis(z)),
42              "large_jump_rate": jump_rate(e,a),
43            })
44        e = EPS
45        def vals(k): return np.array([r[k] for r in rows])
46        obs = {
47          "variance_exponent": slope(e, vals("sigma2_theory")),
48          "naive_error_exponent": slope(e, vals("naive_scale_sigma")),
49          "compensated_W1_proxy_exponent": slope(e, vals("compensated_w1_proxy_rho_over_sigma2")),
50          "residual_normality_exponent": slope(e, vals("normality_proxy_rho_over_sigma3")),
51          "large_jump_rate_exponent": slope(e, vals("large_jump_rate")),
52          "mean_gaussian_relative_variance_error": float(np.mean(vals("variance_relative_error"))),
53        }
54        pred = {
55          "variance_exponent": 2-a,
56          "naive_error_exponent": 1-a/2,
57          "compensated_W1_proxy_exponent": 1.0,
58          "residual_normality_exponent": a/2,
59          "large_jump_rate_exponent": -a,
60        }
61        out["results"][str(a)] = {"observed": obs, "predicted": pred, "rows": rows}
62    with open("results.json", "w") as f: json.dump(out, f, indent=2)
63    for a,d in out["results"].items():
64        print(f"alpha={a}")
65        for k in d["predicted"]:
66            print(f"  {k}: observed={d['observed'][k]:.5f} predicted={d['predicted'][k]:.5f}")
67        print(f"  MC relative variance error={d['observed']['mean_gaussian_relative_variance_error']:.5f}")
68        for r in d["rows"]:
69            print(f"    eps={r['epsilon']:.3g} sigma2={r['sigma2_theory']:.5g} "
70                  f"naive={r['naive_scale_sigma']:.5g} comp_proxy={r['compensated_w1_proxy_rho_over_sigma2']:.5g} "
71                  f"normality={r['normality_proxy_rho_over_sigma3']:.5g} rate={r['large_jump_rate']:.5g}")
72
73if __name__ == "__main__": main()