重标极差（R/S）分析法计算Hurst指数（Python）

记录下自己论文中态势预测问题中，使用重标极差分析法对时间序列数据集进行可预测分析的过程。网上找到的相关R/S计算Hurst指数的代码，大多没有按照标准计算过程来实现，而相关论文中使用Hurst指数时，往往采用了对数散点图的方式来展示。

最终实现效果：

水平有限，只是简单实现，没有采用矩阵等方式加速计算过程。

import numpy as np
import matplotlib.pyplot as plt


def my_hurst(obs: list) -> int:
    """
    严格按照 R/S 标准流程推导
    :param obs: 时间序列数据
    :return: Hurst指数
    """
    x_res = []
    y_res = []
    # 划分时间序列的间隔
    step = 3
    for n in range(step, len(obs), step):
        x_res.append(np.log(n))
        # 按长度n将时间序列分成子序列
        sub_obs = []
        for i in range(0, len(obs), n):
            sub_obs.append(obs[i:i + n])
        # 每个子序列的均值
        sub_mean = []
        # 计算每个子序列的标准差
        sub_S = []
        # 计算每个子序列的极差
        sub_R = []
        for sub in sub_obs:
            # 每个子序列的累计离差
            sub_cumulative_deviation = []
            mean = sum(sub) / n
            sub_mean.append(mean)
            for i in range(n):
                deviation = [x - mean for x in sub[:i]]
                sub_cumulative_deviation.append(sum(deviation))
            sub_R.append(max(sub_cumulative_deviation) - min(sub_cumulative_deviation))
            # 离差平方
            z = [(x - mean) ** 2 for x in sub]
            sub_S.append(np.sqrt((sum(z) + 0.001) / n))

        # 计算每组RS值
        sub_RS = []
        for i in range(len(sub_S)):
            R = sub_R[i]
            S = sub_S[i]
            sub_RS.append(R / S)
        # 计算重标极差
        RS = sum(sub_RS) / len(sub_RS)
        y_res.append(np.log(RS))

    plt.xlabel("ln(n)")
    plt.ylabel("ln(R/S)")
    plt.xlim(1, 5)  # 制定坐标轴范围大小
    plt.ylim(0, 5)

    # 拟合
    h, b = np.polyfit(x_res, y_res, 1)
    x = [0] + x_res + [5]
    x = np.array(x)
    y = h * x + b

    plt.scatter(x_res, y_res, label='Hurst R/S')
    if b > 0:
        lable = "y = %.3f * x + %.3f" % (h, b)
    else:
        lable = "y = %.3f * x - %.3f" % (h, abs(b))
    plt.plot(x, y, label=lable)
    plt.legend(loc='best')
    plt.show()
    return h

运行结果：

Hurst =  0.8170858515778376