时间序列分析之 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个连续子区间，中每一点记作 $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}}$ 。

$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）计算各子区间上的重标极差值（即：对子区间上极差重新标度）：

，

其中， $\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）重复以上过程，并将所得的平均重标极差值对n进行双对数回归：

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

可得到相应的 Hurst指数。

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

参见文章相关链接：西北大学研究团队基于Sentinel-1A卫星影像完成兰新高速铁路沿线部分区域滑坡监测-ESPL

5. matlab 代码

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

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

打开下面的网址查看：file:///d:/Downloads/2000%E2%80%%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..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..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 @.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/

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和负值，得出该值，时序选用其它分析方法。

免责声明：本站所有文章内容,图片，视频等均是来源于用户投稿和互联网及文摘转载整编而成，不代表本站观点，不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益，请在线联系站长,一经查实,本站将立刻删除。本文来自网络,若有侵权，请联系删除，如若转载，请注明出处：https://haidsoft.com/156399.html

时间序列分析之 Hurst 指数

1. Hurst 指数介绍

2. 形式

3. 计算方法

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

5. matlab 代码

6. python 实现

相关推荐

发表回复