示例08 时间依赖 Hurst 指数 HT

"""ex08_ht.py — 第 7 章 时间依赖 Hurst 指数 HT

目标：用 fathon.HT 追踪 Hurst 指数随时间的变化，
构造一个「前半持久、后半反持久」的合成序列，验证 HT 能识别制度切换。

同时演示用预计算的 MFDFA q=0 拟合（q0Fit）加速 HT 计算。

对应教程：tutorials/07-ht.md。

运行：
python ex08_ht.py
"""

import time
import numpy as np
import fathon
from fathon import fathonUtils as fu
import matplotlib.pyplot as plt

def build_regime_switching(L, seed=42):
"""前半白噪声（持久），后半反相关（反持久）。"""
rng = np.random.default_rng(seed)
inc1 = rng.standard_normal(L // 2)
inc2 = -0.5 * inc1 + rng.standard_normal(L // 2) * 0.5
return np.concatenate([inc1, inc2])

def plot_ht(profile, ht, scales, L):
ht_arr = np.asarray(ht)

fig, axes = plt.subplots(2, 1, figsize=(12, 6))

axes[0].plot(profile, lw=0.5, color='steelblue')
axes[0].axvline(L // 2, color='red', ls='--', alpha=0.7,
                label='制度切换点')
axes[0].set_ylabel('Y(t)')
axes[0].set_title('Profile')
axes[0].legend()

if ht_arr.ndim == 2:
    for i, s in enumerate(scales):
        axes[1].plot(ht_arr[i], lw=0.8, label=f's={s}')
else:
    axes[1].plot(ht_arr, lw=0.8, label='H(t)')
axes[1].axhline(0.5, color='gray', ls=':', alpha=0.7)
axes[1].axvline(L // 2, color='red', ls='--', alpha=0.7)
axes[1].set_ylabel('局部 Hurst')
axes[1].set_xlabel('t')
axes[1].set_title('时间依赖 Hurst 指数')
axes[1].legend()
axes[1].set_ylim(0, 1.2)

plt.tight_layout()
​

def main():
np.random.seed(42)
L = 8000
increments = build_regime_switching(L)
profile = fu.toAggregated(increments)

# ── 路径 A：直接 computeHt（内部会跑 MFDFA q=0）──
scales = [100, 200, 500]
ht_obj = fathon.HT(profile)
t0 = time.time()
ht = ht_obj.computeHt(scales, polOrd=1, mfdfaPolOrd=1, verbose=True)
t_direct = time.time() - t0
print(f"直接 computeHt 耗时：{t_direct:.1f}s")

plot_ht(profile, ht, scales, L)

# ── 路径 B：预计算 q0Fit 复用 ──
wins = fu.linRangeByStep(10, L // 4)
q_list = np.arange(-2, 3, 0.5)

mfdfa = fathon.MFDFA(profile)
mfdfa.computeFlucVec(wins, q_list, revSeg=True)
h_q, h_intercept = mfdfa.fitFlucVec(logBase=np.e)
q0_idx = int(np.argmin(np.abs(q_list - 0)))
q0Fit = [h_q[q0_idx], h_intercept[q0_idx]]

ht_obj2 = fathon.HT(profile)
t0 = time.time()
ht2 = ht_obj2.computeHt([200, 500], polOrd=1, q0Fit=q0Fit)
t_reuse = time.time() - t0
print(f"复用 q0Fit 耗时：{t_reuse:.1f}s  (省去内部 MFDFA)")

plt.show()
print("\n观察：H(t) 在切换点附近应出现明显变化（前半偏大、后半偏小）。")
​

if name == "main":
main()