常数/零值时间序列的Hurst系数计算问题及方法验证

问题 #

我需要计算一批（batch）时间序列的Hurst系数，这些序列来自t=0时刻受扰动后自由运行2000步的系统，要为每个系统参数匹配对应的H值。但在处理零值和常数型时间序列时遇到了问题：

我尝试了两种Hurst指数计算方法：

1. 基于维基百科的R/S比率法，代码如下：

import numpy as np

n_min = 10
ns = np.arange(n_min, len(time_series), 5)

R = []
S = []
E = []
for n in ns:
    Xt = time_series[:n]
    m = np.mean(Xt)
    Yt = Xt - m
    Zt = [sum(Yt[:i]) for i in range(len(Yt))]
    # Compute the range
    R.append(max(Zt)-min(Zt))
    # Compute standard deviation
    S.append(np.std(Xt))
    E.append(np.mean(np.array(R)/np.array(S)))   

reg = np.polyfit(np.log(ns), np.log(E), 1)
H = reg[0]

2. 基于Stack Overflow的一阶指数法，代码如下：

from numpy import subtract, sqrt, std, polyfit, log10

def hurst(p):  
    tau = []; lagvec = []  
    #  Step through the different lags  
    for lag in range(2,20):  
        #  produce price difference with lag  
        pp = subtract(p[lag:],p[:-lag])  
        #  Write the different lags into a vector  
        lagvec.append(lag)  
        #  Calculate the variance of the differnce vector  
        tau.append(sqrt(std(pp)))  
    #  linear fit to double-log graph (gives power)  
    m = polyfit(log10(lagvec),log10(tau),1)  
    # calculate hurst  
    hurst = m[0]*2  
    # plot lag vs variance  
    #py.plot(lagvec,tau,'o'); show()  
    return hurst  

两种方法得到的结果差异极大：

Zero time series
H (S/R ratio) = 0.4732894
H (1st exponent) = 0.001971641

注：常数时间序列的结果与零值序列几乎一致

由此我对整体时间序列分析的准确性产生疑问，特提出以下问题：

1. 已知Hurst系数适用于平稳分布，能否计算非平稳序列的H值？
2. 若限定分析平稳序列，如何正确计算常数序列的H值？
3. 如何用Python验证结果的统计显著性，以判断哪种实现正确？

回答 #

1. 非平稳序列的Hurst系数计算 #

Hurst系数并非只能用于平稳序列，经典R/S分析最初就是为识别长期记忆性（包括非平稳的分形布朗运动）设计的，但需注意：

• 传统R/S对确定性趋势敏感，会高估H值，需先做去趋势处理（如去趋势R/S分析，DFA）；
• 一阶指数法通过分析序列自相似性计算H，差分操作可消除部分趋势，但面对多项式型趋势，需调整差分阶数或先去趋势。

2. 常数序列的Hurst系数计算 #

常数/零序列是极端平稳且无波动的情况，理论上Hurst系数为1（完全可预测，长期记忆性最强），但常规计算会遇到除以0的问题，需特殊处理：

• R/S方法：先判断序列是否为常数（np.all(time_series == time_series[0])），若是直接返回H=1；否则按正常流程计算，同时跳过S=0的情况；
• 一阶指数法：同样先判断常数序列，直接返回H=1；若序列因浮点误差接近常数，需跳过标准差为0的差分结果，避免log(0)报错。

另外，你的第一种方法代码存在逻辑错误：R/S分析需对每个窗口n计算单独的R(n)/S(n)，再对log(n)和log(R(n)/S(n))做线性回归，而非先求R/S的均值再回归，这是导致常数序列结果异常的核心原因。修正后的R/S代码如下：

import numpy as np

def hurst_rs_correct(time_series):
    n_min = 10
    ns = np.arange(n_min, len(time_series), 5)
    rs_vals = []
    valid_ns = []
    for n in ns:
        Xt = time_series[:n]
        m = np.mean(Xt)
        Yt = Xt - m
        Zt = np.cumsum(Yt)  # 用numpy内置函数替代循环，更高效
        R_val = max(Zt) - min(Zt)
        S_val = np.std(Xt, ddof=1)  # 使用样本标准差，更符合统计定义
        if S_val == 0:
            continue
        rs_vals.append(R_val / S_val)
        valid_ns.append(n)
    if len(rs_vals) < 2:
        return np.nan  # 有效窗口不足，返回NaN
    log_n = np.log(valid_ns)
    log_rs = np.log(rs_vals)
    reg = np.polyfit(log_n, log_rs, 1)
    return reg[0]

3. Python验证结果统计显著性的方法 #

（1）蒙特卡洛模拟对比已知H值的序列

生成具有确定Hurst系数的模拟序列（如分形布朗运动），用两种方法计算并对比误差：

from fbm import FBM  # 需安装：pip install fbm

# 生成已知H=0.7的分形布朗运动序列
fbm_gen = FBM(n=2000, hurst=0.7, length=1, method='daviesharte')
fbm_series = fbm_gen.fbm()

# 调用修正后的R/S方法和一阶指数法
h_rs = hurst_rs_correct(fbm_series)
h_fo = hurst(fbm_series)

print(f"真实H值：0.7，修正后R/S方法H值：{h_rs:.4f}，一阶指数法H值：{h_fo:.4f}")

多次运行模拟，统计两种方法的均值误差、标准差，判断哪种更接近真实值。

（2）常数序列的显著性验证

• 给常数序列添加微小高斯噪声，生成近常数序列，观察两种方法的H值是否趋近于理论值1；
• 置换检验：打乱原序列（破坏相关性），计算打乱后的H值，对比原序列H值是否显著不同，判断方法的有效性。

（3）代码逻辑验证

检查两种方法的计算流程是否符合理论定义：

• R/S方法需确保每个窗口n对应单独的R/S值，回归对象是log(n)与log(R/S)的成对数据；
• 一阶指数法需确保差分后的标准差不为0，避免无效计算。

内容的提问来源于stack exchange，提问作者Emmanuel Calvet