时间序列分析之 Hurst 指数

顺心而行...

已于 2024-05-14 11:14:23 修改

阅读量2k

点赞数 6

文章标签：算法经验分享笔记

于 2024-05-10 11:54:28 首次发布

本文链接：https://blog.csdn.net/weixin_41652574/article/details/138655743

版权

1. Hurst 指数介绍

Hurst指数，又称为赫斯特指数或赫斯特维特指数，是一种用于分析时间序列数据的统计指标。

它最早由机械工程师Harold Edwin Hurst于1951年引入，并用于研究尼罗河洪水的周期性。

Hurst指数可以帮助我们理解时间序列数据的长期趋势，以及预测未来趋势的稳定性。

2. 形式

赫斯特指数有三种形式：

1．如果H=0.5，表明时间序列接近无定向运动的随机游走；

2．如果0.5<H<1，表明时间序列存在长期记忆性，未来的趋势和过去的趋势相关，继续保持现有趋势的可能性强；

3．如果0≤H<0.5，表明粉红噪声(反持续性)即均值回复过程。即很有可能是记忆的转弱，趋势结束和反转的开始。

3. 计算方法

HURST指数的计算方法主要有七种：

聚合方差法(Aggregated Variance method)，R/S分析法(R/S method)，周期图法(Periodogram method)，绝对值法(Absolute Value method)，残差方差法(Variance of residuals)，小波分析法(Abry-Veitch method)，Whittle法(Whittle estimator)。

R/S分析法，即重标极差分析法。用此法计算HURST指数，不仅计算量大，且方法繁杂。一般都是针对少数代表性指数，且多半是用月（周）数据分析的。（可用于形变序列分析）

Hurst指数的计算方式相对简单。以最常用的重标极差法为例，其计算过程如下：

对于Hurst指数的计算方法(重标极差法)，其计算流程一般分为三步：
1、将长度为N的时间序列划分为长度为n的A个连续子区间 $I_a(a=1,......,A)$ ， $I_a$ 中每一点记作 $R_{k,a}$ ；

2、针对不同的子区间长度n，计算A个子区间的平均重标极差：

$\left(\mathrm{R/S}\right)_{\mathrm{n}}=\frac{1}{\mathrm{A}}\sum_{\mathrm{k=1}}^{\mathrm{A}}(\mathrm{R/S})_{\mathrm{k}}$ 。

其中每个子区间上重标极差的计算步骤：
1）对每个长度为n的子区间，在对相应时间序列作零均值化处理后，计算其累计离差：

$X_{N,a}=\sum_{\mathrm{k=1}}^{\mathrm{N}}(\mathrm{R_{k,a}}-\frac{1}{\mathrm{N}}\sum_{\mathrm{k=1}}^{\mathrm{n}}\mathrm{R_{k,a}})$

2）定义单个子区间上的极差：

$\mathrm{R_{a}=max\bigl(X_{k,a}\bigr)-min\bigl(X_{k,a}\bigr),k=1,2,......,n}$

3）计算各子区间上的重标极差值（即：对子区间上极差重新标度）：

$(R/S)_a=R_a/S_a$ ，

其中， $\mathrm{S_{i}=\sqrt{\frac{1}{n}\sum_{k=1}^{n}(R_{k,i}-\frac{1}{N}\sum_{k=1}^{n}R_{k,a})^{2}} , i=1,2,......,A}$

3、由于样本的平均重标极差值与样本长度之间存在标度关系，即：

$\mathrm{(R/S)}_{\mathrm{n}}=\mathrm{c\times n^{H}}$

因此，对不同时间尺度（即：不同划分长度n）重复以上过程，并将所得的平均重标极差值 $(R/S)_n$ 对n进行双对数回归：

$\mathrm{log(R/S)_n=log(C)+H\cdot log(n)}$

可得到相应的 Hurst指数。

4. Hurst 指数在形变分析中的应用

5. matlab 代码

【B站有视频：https://www.bilibili.com/video/BV19C411L7K4/?】

视频中的参考文献：《2000—2020 年黑龙江省植被时空变化对气候因子响应》

打开下面的网址查看：file:///d:/Downloads/2000%E2%80%942020%E5%B9%B4%E9%BB%91%E9%BE%99%E6%B1%9F...%E6%A4%8D%E8%A2%AB%E6%97%B6%E7%A9%BA%E5%8F%98%E5%8C%96%E5%AF%B9%E6%B0%94%E5%80%99%E5%9B%A0%E5%AD%90%E5%93%8D%E5%BA%94_%E5%88%98%E6%99%BA%E6%BA%90.pdf

（1）多年NDVI数据的HURST指数计算

参考：https://mp.weixin.qq.com/s/pXXwRU9XMdom10HoLx45uA

clear
[aa,R]=geotiffread('H:\Global\NDVI3g\GIMMSraster\raster\最大合成\GIMMS_NDVI2015.tif');%先投影信息
info=geotiffinfo('H:\Global\NDVI3g\GIMMSraster\raster\最大合成\GIMMS_NDVI2015.tif');
ndvisum=zeros(size(aa,1)*size(aa,2),34);%34期数据
for year=1982:2015
    filename=strcat('H:\Global\NDVI3g\GIMMSraster\raster\最大合成\GIMMS_NDVI',int2str(year),'.tif');
    ndvi=importdata(filename);
    ndvi=reshape(ndvi,size(ndvi,1)*size(ndvi,2),1);
    ndvisum(:,year-1981)=ndvi;
end
hsum=zeros(size(aa,1),size(aa,2))+NaN;
for kk=1:size(ndvisum,1);
    ndvi=ndvisum(kk,:);
    if min(ndvi)>0
        ndvi_cf=[];
        for i=1:length(ndvi)-1
            ndvi_cf1=ndvi(i+1)-ndvi(i);
            ndvi_cf=[ndvi_cf,ndvi_cf1];
        end
        M=[];
        for i=1:size(ndvi_cf,2)
            M1=mean(ndvi_cf(1:i));
            M=[M,M1];
        end
        S=[];

        for i=1:size(ndvi_cf,2)
            S1=std(ndvi_cf(1:i))*sqrt((i-1)/i);
            S=[S,S1];
        end

        for i=1:size(ndvi_cf,2)
            for j=1:i
                der(j)=ndvi_cf(1,j)-M(1,i);
                cum=cumsum(der);
                RR(i)=max(cum)-min(cum);
            end
        end

        RS=S(2:size(ndvi_cf,2)).\RR(2:size(ndvi_cf,2));
        T=[];
        for i=1:size(ndvi_cf,2)
            T1=i;
            T=[T,T1];
        end
        lag=T(2:size(ndvi_cf,2));                  
        g=polyfit(log(lag/2),log(RS),1);                
        H=g(1);
        hsum(kk)=H;
        clear der
    end
end
geotiffwrite('1982-2015年全球NDVI_Hurst指数.tif',hsum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

（2）多年NDVI数据的sen趋势分析

参考：https://mp.weixin.qq.com/s/Rc4j-FRP41_TWozczemrHA

Theil-Sen median趋势分析是一种稳健的非参数计算方法，计算式如下：

$S_{ET}=median\left(\frac{ET_j-ET_i}{j-i}\right), 2000\leqslant i<j\leqslant2014$

公式中 $S_{ET}$ 指计算n(n-1)/2 个数据组合的斜率的中位数， $ET_j$ 和 $ET_j$ 代表i和j年的ET值。如果 $S_{ET}$ ＞0，则ET呈上升趋势，否则，ET呈下降趋势。

Matlab具体代码如下：

% @author 1154318421@qq.com
[a,R]=geotiffread('E:\GIS\NDVI\2000.tif');%先导入投影信息
info=geotiffinfo('E:\GIS\NDVI\2000.tif');
[m,n]=size(a);
cd=2020-2000+1;%时间跨度，根据需要自行修改
datasum=zeros(m*n,cd)+NaN; 
k=1;
for year=2000:2020 %起始年份
    filename=['E:\GIS\ENVI\xinjiang',int2str(year),'2000.tif'];
    data=importdata(filename);
    data=reshape(data,m*n,1);
    datasum(:,k)=data;
    k=k+1;
end
result=zeros(m,n)+NaN;
for i=1:size(datasum,1)
    data=datasum(i,:);
    if min(data)>0 %判断是否是有效值,我这里的有效值必须大于0
        valuesum=[];
        for k1=2:cd
            for k2=1:(k1-1)
                cz=data(k1)-data(k2);
                jl=k1-k2;
                value=cz./jl;
                valuesum=[valuesum;value];
            end
        end
        value=median(valuesum);
        result(i)=value;
    end
end
filename=['E:\GIS\基于SEN的NDVI变化趋势.tif'];
geotiffwrite(filename,result,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)

6. python 实现

（1）在python中，可以使用 pyeeg 库实现 hurst 指数计算

安装 pyeeg 库：

pip install pyeeg -i https://pypi.tuna.tsinghua.edu.cn/simple/ --trusted-host tuna.tsinghua.edu.cn

导入与调用：

import pyeeg

hurst = pyeeg.hurst(data) # data 为时间序列

实现过程参考：https://www.cnblogs.com/ningjing213/p/10745772.html

def Hurst(data):
    n = 6
    data = pd.Series(data).pct_change()[1:]
    ARS = list()
    lag = list()
    for i in range(n):
        m = 2 ** i
        size = np.size(data) // m
        lag.append(size)
        panel = {}
        for j in range(m):
            panel[str(j)] = data[j*size:(j+1)*size].values

        panel = pd.DataFrame(panel)
        mean = panel.mean()
        Deviation = (panel - mean).cumsum()
        maxi = Deviation.max()
        mini = Deviation.min()
        sigma = panel.std()
        RS = maxi - mini
        RS = RS / sigma
        ARS.append(RS.mean())

    lag = np.log10(lag)
    ARS = np.log10(ARS)
    hurst_exponent = np.polyfit(lag, ARS, 1)
    hurst = hurst_exponent[0]

    return hurst

（2）其它参考：https://cloud.tencent.com/developer/article/1518357

def hurst(ts, if_detail=False):
    N = len(ts)
    if N < 20:
        raise ValueError("Time series is too short! input series ought to have at least 20 samples!")
    max_k = int(np.floor(N / 2))
    n_all = []
    RS_all = []
    # k 是子区间长度
    for k in range(5, max_k + 1):
        subset_list = np.array(ts[:N - N % k]).reshape(-1, k).T
        # 累积极差
        cumsum_list = (subset_list - subset_list.mean(axis=0)).cumsum(axis=0)
        R = cumsum_list.max(axis=0) - cumsum_list.min(axis=0)
        S = (((subset_list - subset_list.mean(axis=0)) ** 2).mean(axis=0)) ** 0.5
        RS = (R / S).mean()
        n_all.append(k)
        RS_all.append(RS)
    # print(k)
    # R_S_all = pd.DataFrame(R_S_all)
    # R_S_all['logN'] = np.log10(R_S_all.n)
    # R_S_all['logRS'] = np.log10(R_S_all.RS)
    Hurst_exponent = np.polyfit(np.log10(n_all), np.log10(RS_all), 1)[0]
    if if_detail:
        n_all = np.array(n_all)
        RS_all = np.array(RS_all)
        Vn = RS_all / np.sqrt(n_all)
        res = pd.DataFrame([n_all, RS_all, Vn]).T
        res.columns = ['n', 'RS', 'Vn']
        return res, Hurst_exponent
    # plt.plot(R_S_all.logN,R_S_all.logRS)
    else:
        return Hurst_exponent

（3）上述python代码运行都存在超过 1 的 hurst 值，有人说属于正常现象，超过1的话序列非平穏。以下是另一参考，避免了 hurst 大于1的情况。

参考：https://en.wikipedia.org/wiki/Hurst_exponent

 # -*- coding: utf-8 -*-
    # Reference: https://en.wikipedia.org/wiki/Hurst_exponent
    # python 3.6.2 AMD64
    # 2018/4/19
    # Calculate Hurst exponent based on Rescaled range (R/S) analysis
    # How to use (example):
    # import Hurst 
    # ts = list(range(50))
    # hurst = Hurst.hurst(ts)
    # Tip: ts has to be object list(n_samples,) or np.array(n_samples,)

    __Author__ = "Ryan Wang"


    def hurst(ts):
        ts = list(ts)
        N = len(ts)
        if N < 20:
            raise ValueError("Time series is too short! input series ought to have at least 20 samples!")

        max_k = int(np.floor(N / 2))
        R_S_dict = []
        for k in range(10, max_k + 1):
            R, S = 0, 0
            # split ts into subsets
            subset_list = [ts[i:i + k] for i in range(0, N, k)]
            if np.mod(N, k) > 0:
                subset_list.pop()
                # tail = subset_list.pop()
                # subset_list[-1].extend(tail)
            # calc mean of every subset
            mean_list = [np.mean(x) for x in subset_list]
            for i in range(len(subset_list)):
                cumsum_list = pd.Series(subset_list[i] - mean_list[i]).cumsum()
                R += max(cumsum_list) - min(cumsum_list)
                S += np.std(subset_list[i])
            R_S_dict.append({"R": R / len(subset_list), "S": S / len(subset_list), "n": k})
        log_R_S = []
        log_n = []
        # print(R_S_dict)
        for i in range(len(R_S_dict)):
            R_S = (R_S_dict[i]["R"] + np.spacing(1)) / (R_S_dict[i]["S"] + np.spacing(1))
            log_R_S.append(np.log10(R_S))
            log_n.append(np.log10(R_S_dict[i]["n"]))
        Hurst_exponent = np.polyfit(log_n, log_R_S, 1)[0]
        return Hurst_exponent

上述代码也并没有解决 hurst 值超过 1 且有负值的情况 ...

查找原因可能是由于数据特性导致：某些类型的时间序列数据可能具有非典型的统计特性，导致Hurst指数超出了0到1的范围。例如，具有非线性、非平稳或非高斯分布的数据可能会导致计算结果偏离预期范围。在这种情况下，可能需要进一步研究数据的性质，并考虑使用其他方法进行分析。

建议：Hurst 计算代码选用（2），不必纠结超过1和负值，得出该值，时序选用其它分析方法。