07 时间依赖 Hurst 指数 HT

第 7 章 时间依赖 Hurst 指数 HT

全局 DFA 给出一个平均 Hurst 指数。当序列的持久性随时间或尺度变化时，HT（Time-dependent Hurst exponent） 提供局部、时变的标度刻画。fathon 实现基于 Ihlen (2012) 的方法。

7.1 何时需要 HT？

HT 在重叠滑动窗口内计算局部波动，并结合 MFDFA 在 q→0 处的标度，得到每个时间点、每个尺度上的局部 Hurst。

7.2 Ihlen 方法直觉

对尺度 s（窗口长度）：

1. 在时刻 t 附近取长度为 s 的窗口
2. 窗口内对 profile 做多项式去趋势
3. 计算局部波动 F²(t, s)
4. 对给定 s，在 t 上滑动得到 F² 序列
5. 结合 MFDFA q=0 的标度关系，转换为局部 Hurst H(t, s)

fathon 将上述流程封装为 HT.computeHt。

7.3 fathon.HT API

构造

import fathon

ht_obj = fathon.HT(tsVec)  # profile，同 DFA 输入
​

computeHt

ht = ht_obj.computeHt(
    scales,          # 窗口尺度列表，如 [100, 200, 1000]
    polOrd=1,        # 局部去趋势阶数
    mfdfaPolOrd=1,   # MFDFA q=0 拟合的多项式阶数
    q0Fit=[],        # 可选：预计算的 [hq0, hq0_intercept]
    verbose=False
)
​

返回值 ht：时间依赖 Hurst 数组（具体 shape 取决于 scales 数量与序列长度，见下方示例）。

7.4 完整可运行示例

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

np.random.seed(42)
L = 8000

# 构造「前半持久、后半反持久」的合成序列
inc1 = np.random.randn(L // 2)
inc2 = -0.5 * inc1 + np.random.randn(L // 2) * 0.5  # 反相关于前半
increments = np.concatenate([inc1, inc2])
profile = fu.toAggregated(increments)

scales = [100, 200, 500]

ht_obj = fathon.HT(profile)
ht = ht_obj.computeHt(scales, polOrd=1, mfdfaPolOrd=1, verbose=True)

print("HT 结果类型:", type(ht), "shape:", np.asarray(ht).shape)

# ── 可视化 profile 与 HT ──
fig, axes = plt.subplots(2, 1, figsize=(12, 6), sharex=False)

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()

# ht 可能是 (len(scales), n_points) 或 (n_points,) — 按实际 shape 适配
ht_arr = np.asarray(ht)
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()
plt.show()
​

7.5 预计算 q0Fit 加速

若已对同一 profile 跑过 MFDFA，可提取 q=0 的 Hurst 与截距传入，避免 HT 内部重复 MFDFA：

import numpy as np
import fathon
from fathon import fathonUtils as fu

profile = fu.toAggregated(np.random.randn(6000))
wins = fu.linRangeByStep(10, 1500)
q_list = np.arange(-2, 3, 0.5)

# 先 MFDFA
mfdfa = fathon.MFDFA(profile)
n, F = mfdfa.computeFlucVec(wins, q_list, revSeg=True)
h_q, h_intercept = mfdfa.fitFlucVec(logBase=np.e)

q0_idx = np.argmin(np.abs(q_list - 0))
q0Fit = [h_q[q0_idx], h_intercept[q0_idx]]

# HT 复用 q=0 拟合
ht_obj = fathon.HT(profile)
ht = ht_obj.computeHt([200, 500], polOrd=1, q0Fit=q0Fit)
print("HT 计算完成，使用预计算 q0Fit")
​

注意：q0Fit 中的截距必须来自自然对数（logBase=e） 的拟合，与 fitFlucVec(logBase=np.e) 一致。

7.6 HT 结果解读

7.7 HT vs 滑动窗口 DFA

7.8 计算成本与尺度选择

HT 比单次 DFA 显著更慢（多尺度 × 滑动 × 内部 MFDFA）。建议：

• 序列长度 L ≥ 5000 再考虑 HT
• scales 选 2～4 个代表性尺度，勿过多
• 先用全局 DFA/MFDFA 确认存在标度，再跑 HT
• 用 q0Fit 避免重复 MFDFA

7.9 saveObject

与 DFA 相同：

ht_obj.saveObject("ht_result")
​

7.10 本章小结

profile ──► HT ──► computeHt(scales) ──► H(t, s)
                      │
              可选 q0Fit ◄── MFDFA q=0
​

HT 适合非平稳 persistence 与制度切换检测；平稳序列用全局 DFA 即可。

7.11 动手实验

1. 对纯白噪声 profile 跑 HT，观察 H(t) 是否在 0.5 附近波动。
2. 比较 scales=[50] 与 scales=[50,200,500] 的结果差异。
3. 用第 7.5 节流程，对比有无 q0Fit 的运行时间（verbose=True）。

7.12 配套可运行示例

「时间依赖 Hurst（HT）」示例完整实现本章流程：构造前半持久、后半反持久的制度切换序列，对比「直接 computeHt」与「预计算 q0Fit 复用」两条路径的耗时。下面是可直接运行的完整脚本：

"""时间依赖 Hurst 指数（HT）完整示例。

构造一个「前半持久、后半反持久」的合成序列，用 HT 追踪 Hurst 随时间
的变化，验证它能识别制度切换；并演示用预计算的 q0Fit 加速。
依赖：pip install fathon numpy
"""
import time
import numpy as np
import fathon
from fathon import fathonUtils as fu


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 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=False)
    t_direct = time.time() - t0
    ht_arr = np.asarray(ht)
    print(f"直接 computeHt 耗时：{t_direct:.1f}s")
    print(f"HT 结果 shape：{ht_arr.shape}")

    # ── 路径 B：预计算 q0Fit 复用，省去内部 MFDFA ──
    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)")

    # ── 解读：观察 H(t) 在切换点附近的变化 ──
    if ht_arr.ndim == 2:
        for i, s in enumerate(scales):
            half = ht_arr[i].shape[0] // 2
            h_before = np.nanmean(ht_arr[i][:half])
            h_after = np.nanmean(ht_arr[i][half:])
            print(f"  scale={s}: 切换前 H≈{h_before:.3f}  切换后 H≈{h_after:.3f}")
    print("切换点两侧 H 应有明显变化——这是全局 DFA 看不到的局部信息。")


if __name__ == "__main__":
    main()
​

💡 注意 q0Fit 中的截距必须来自自然对数（logBase=e）的拟合，与 fitFlucVec(logBase=np.e) 一致，否则 HT 结果会有偏差。

下一章集中讲解 fathonUtils 与窗口设计的工程细节。