1"""Tempered-stable volatility clock MVP and numerical verification."""
  2import json, math
  3from pathlib import Path
  4import numpy as np
  5
  6
  7def positive_stable_laplace_one(alpha, rng, size):
  8    """Kanter sampler: positive stable S with E exp(-u S)=exp(-u**alpha)."""
  9    u = rng.uniform(1e-9, math.pi - 1e-9, size)
 10    e = rng.exponential(1.0, size)
 11    return (np.sin(alpha*u) / np.sin(u)**(1.0/alpha) *
 12            (np.sin((1-alpha)*u) / e)**((1-alpha)/alpha))
 13
 14
 15def ts_sample(alpha, theta, delta, rng, size):
 16    """Exact rejection sampler for TS(alpha,theta,delta).
 17
 18    If S has Laplace exp(-s^alpha), delta^(1/alpha)S exponentially tilted
 19    by exp(-theta*x) has the requested tempered-stable law.
 20    """
 21    scale = delta ** (1.0/alpha)
 22    out = np.empty(size, dtype=float)
 23    filled = 0
 24    while filled < size:
 25        # Acceptance is exp(-delta*theta**alpha), generally high for settings here.
 26        n = max(128, int((size-filled) / max(math.exp(-delta*theta**alpha), .05) * 1.15))
 27        x = scale * positive_stable_laplace_one(alpha, rng, n)
 28        keep = rng.random(n) < np.exp(-theta*x)
 29        got = x[keep]
 30        take = min(got.size, size-filled)
 31        if take:
 32            out[filled:filled+take] = got[:take]
 33            filled += take
 34    return out
 35
 36
 37def params(alpha, theta, phi):
 38    delta = (1-phi) * theta**(1-alpha) / alpha
 39    v = (1-alpha)/(theta*(1+phi))
 40    return delta, v
 41
 42
 43def simulate(alpha, theta, phi, n_paths, length, rng, burn=250):
 44    delta, _ = params(alpha, theta, phi)
 45    a = np.ones(n_paths)
 46    # Burn-in makes the finite recursion close to stationary; mean initialization
 47    # also avoids an unnecessarily long burn for phi near one.
 48    for _ in range(burn):
 49        a = phi*a + ts_sample(alpha, theta, delta, rng, n_paths)
 50    A = np.empty((n_paths, length))
 51    for i in range(length):
 52        a = phi*a + ts_sample(alpha, theta, delta, rng, n_paths)
 53        A[:, i] = a
 54    z = rng.standard_normal((n_paths, length))
 55    return A, np.sqrt(A)*z
 56
 57
 58def excess_kurtosis(x):
 59    x = x.ravel()
 60    m2 = np.mean(x*x)
 61    return np.mean(x**4)/(m2*m2)-3
 62
 63
 64def lag_sq_acf(x, lag):
 65    q = x*x
 66    q = q - q.mean()
 67    return np.mean(q[:, :-lag]*q[:, lag:]) / np.mean(q*q)
 68
 69
 70def laplace_check(alpha, theta, delta, rng, n=200000):
 71    x = ts_sample(alpha, theta, delta, rng, n)
 72    rows=[]
 73    for u in (0.2, 1.0, 3.0):
 74        empirical=np.mean(np.exp(-u*x))
 75        theory=math.exp(-delta*((theta+u)**alpha-theta**alpha))
 76        rows.append({'u':u, 'empirical':empirical, 'theory':theory,
 77                     'abs_error':abs(empirical-theory)})
 78    return rows
 79
 80
 81def run():
 82    rng=np.random.default_rng(1142)
 83    # Prediction 1: stationary mean is one across persistence values.
 84    mean_sweep=[]
 85    # Prediction 2: variance / excess kurtosis follows closed form across alpha.
 86    var_sweep=[]
 87    # Prediction 3: squared-return ACF equals phi^h*v/(2+3v), with geometric decay.
 88    acf_sweep=[]
 89    for phi in (0.0, 0.3, 0.7, 0.9):
 90        alpha,theta=.65,.8
 91        A,r=simulate(alpha,theta,phi,3500,160,rng)
 92        delta,v=params(alpha,theta,phi)
 93        mean_sweep.append({'phi':phi,'pred_mean':1.0,'obs_mean':float(A.mean()),
 94                           'pred_var':v,'obs_var':float(A.var())})
 95        acf_sweep.append({'phi':phi,'pred_lag1':phi*v/(2+3*v),
 96                          'obs_lag1':float(lag_sq_acf(r,1)),
 97                          'pred_lag5':phi**5*v/(2+3*v),
 98                          'obs_lag5':float(lag_sq_acf(r,5))})
 99    for alpha in (.35,.55,.75):
100        theta=.8; phi=.7
101        A,r=simulate(alpha,theta,phi,3500,160,rng)
102        delta,v=params(alpha,theta,phi)
103        var_sweep.append({'alpha':alpha,'pred_var':v,'obs_var':float(A.var()),
104                          'pred_K':3*v,'obs_K':float(excess_kurtosis(r)),
105                          'delta':delta})
106    # Mechanism comparison at identical dimensions: iid Gaussian has no clustering.
107    alpha,theta,phi=.65,.8,.7
108    A,r=simulate(alpha,theta,phi,3500,160,rng)
109    iid=rng.standard_normal(r.shape)
110    delta,v=params(alpha,theta,phi)
111    baseline={'K':float(excess_kurtosis(iid)), 'rho1':float(lag_sq_acf(iid,1))}
112    idea={'K':float(excess_kurtosis(r)), 'rho1':float(lag_sq_acf(r,1)),
113          'pred_K':3*v, 'pred_rho1':phi*v/(2+3*v)}
114    result={'laplace':laplace_check(.65,.8,delta,rng),
115            'mean_variance_sweep':mean_sweep,'kurtosis_sweep':var_sweep,
116            'acf_sweep':acf_sweep,'baseline_iid_gaussian':baseline,
117            'clock':idea}
118    Path('results.json').write_text(json.dumps(result,indent=2))
119    print(json.dumps(result,indent=2))
120
121if __name__=='__main__': run()