200字范文 > 采用R/S分析法的Hurst指数估计算法——Python实现

采用R/S分析法的Hurst指数估计算法——Python实现

时间：2019-01-26 03:41:28

文章目录

一、前言1. R/S分析法起源2. Hurst指数定义3. R/S 分析法 Hurst 指数的估计二、算法伪码三、Python代码四、代码测试1. 数据2. 结论五、总结

一、前言

代码附在文末了

1. R/S分析法起源

“Hurst 指数”或“Hurst 系数”由研究员Harold Edwin Hurst在研究罗河旱涝更替的情况时，为研究水利的实际问题发明，以衡量时间序列的长期记忆能力。Hurst 指数又被称为“指数依赖性”或“指数长期依赖性”，它能够量化时间序列的相对趋势。

现在有很多 Hurst 指数估计值的算法，最有名的就是 Mandelbrot 和 Wallis 基于 Hurst 的水利研究结果使用的重标极差 R/S 方法。

2. Hurst指数定义

设有一个序列 X={X1,X2,…}X = \{X_1, X_2, … \}X={X1,X2,…} ，nnn 是观测到的时间跨度，R(n)R(n)R(n) 是前 nnn 个值的取值范围，S(n)S(n)S(n) 是它们的标准差，EEE 符号表示求期望值，CCC 是一个常数。则序列 XXX 的 Hurst指数（后面以 HHH 表示）的原始定义如下式所示：

E(R(n)S(n))=CnH,whilen→∞(1)E(\dfrac{R(n)}{S(n)})\ =\ Cn^H,\ \ \ \ while\ n \rightarrow \infty \tag{1} E(S(n)R(n))=CnH,whilen→∞(1)

3. R/S 分析法 Hurst 指数的估计

利用 R/S 分析法，把一个长度总共为 NNN 的时间序列分成长度分别为 n={N,N2,N4,…}n = \{N, \dfrac{N}{2}, \dfrac{N}{4}, … \}n={N,2N,4N,…} 的短序列。对于每一个 nnn，便可以计算其重标极差 R(n)S(n)\dfrac{R(n)}{S(n)}S(n)R(n)。

又因为数据的 Hurst 指数满足幂定律(1)式，便可以画出 log[R(n)S(n)]log[\dfrac{R(n)}{S(n)}]log[S(n)R(n)] 关于 lognlog\ nlogn 的图形，通过拟合直线的斜率得到 HHH 的值。

二、算法伪码

任意长度序列(n≥8n \ge 8n≥8) Hurst 指数计算流程：

三、Python代码

代码中使用到了numpy库：

pip install -i /pypi/simple numpy

计算 Hurst 指数代码如下：

def Hurst(ts):'''Parameters----------ts : Iterable Object.A time series or a list.Raises------ValueErrorIf input ts is not iterable then raise error.Returns-------H : FloatThe Hurst-index of this series.'''if not isinstance(ts, Iterable):raise ValueError("This sequence is not iterable !")ts = np.array(ts)# N is use for storge the length sequenceN, RS, n = [], [], len(ts)while (True):N.append(n)# Calculate the average value of the seriesm = np.mean(ts)# Construct mean adjustment sequencemean_adj = ts - m# Construct cumulative deviation sequencecumulative_dvi = np.cumsum(mean_adj)# Calculate sequence rangesrange = max(cumulative_dvi) - min(cumulative_dvi)# Calculate the unbiased standard deviation of this sequenceunbiased_std_dvi = np.std(ts)# Calculate the rescaled range of this sequence under n lengthRS.append(srange / unbiased_std_dvi)# While n < 4 then breakif n < 4:break# Rebuild this sequence by half lengthts, n = HalfSeries(ts, n)# Get Hurst-index by fit log(RS)~log(n)H = np.polyfit(np.log10(N), np.log10(RS), 1)[0]return H

其中重新生成序列的代码如下：

def HalfSeries(s, n):'''if length(X) is odd:X <- {(X1 + X2) / 2, ..., (Xn-2 + Xn-1) / 2, Xn}n <- (n - 1) / 2else:X <- {(X1 + X2) / 2, ..., (Xn-1 + Xn) / 2}n <- n / 2return X, n'''X = []for i in range(0, len(s) - 1, 2):X.append((s[i] + s[i + 1]) / 2)# if length(s) is oddif len(s) % 2 != 0:X.append(s[-1])n = (n - 1) // 2else:n = n // 2return [np.array(X), n]