在时间序列分析中，有时候会用到hurst指数，今天分享Hurst指数的计算方法。

1 介绍

本节介绍出自《地理数学方法：基础和应用》一书【第 21 章时间序列的 R/S 分析】

R/S 分析是一种基于长程相关思想的时间序列分析方法。这种方法由 H. E. Hurst 于 1965 年最先提出，后来伴随着非线性理论的发展而成长起来。 Hurst 原本是剑桥大学物理学博士，对埃及尼罗河（ Nile）进行了长达 60 年的观测，记录了尼罗河水位变化及其相关的水文、气候数据，发现了所谓的“ Hurst 现象”，开创了 R/S 时间序列分析技术。因此，有学者指出，Hurst 虽然是物理学出身，但却成了一个名副其实的地理学家。
      其实，埃及人保持了几千年的尼罗河水位记录：某些年份洪水泛滥，另一些年份则水势消减。分形理论的创始人 B.B. Mandelbrot 将这种变化划分为两种效应：诺亚（ Noach）效应和约瑟夫（ Joseph）效应。诺亚和约瑟夫这两个名词都本于《圣经》中的人名。大洪水和诺亚方舟的传说众所周知，约瑟夫的故事在东方影响较小。
      诺亚效应意味着不连续性（ discontinuity）：当一个变量变化时，它几乎可以改变得任意快速。换言之，诺亚效应暗示着间断性——许多变化出人意料，不可能进行连续的曲线分析。
      约瑟夫效应意味着持久性（ persistence）： “ 埃及土地上有 7 年大洪水，此后又将有 7 年干旱。 ”（《圣经》）持续性意味着，一个地方受灾越久，它就越可能继续受灾。早年的俄国农民也曾注意到如下地理规律：干旱的发生一般至少持续三年，丰收的出现也通常是连续三年左右，而干旱和丰收年往往以分段连续的方式交替出现。持久性暗示着趋势性，但不等同于周期性。
      诺亚效应与约瑟夫效应叠加的结果是：自然界的趋势是现实的，不过它们可以同样快的来临和同样快地消失。

所谓 R/S 分析，实际上就是重新标度的极差分析（ rescaled range nalysis），简称重标极差分析。具体过程如下。考虑一个时间序列{ξ (t) }，这里ξ (t) = B(t) − B(t −1) ， B(t)为时刻t 的观测值（ t=1,2, …）。对于任意正整数 τ，定义均值序列：

这里τ =1,2,… 代表时滞。用 X(t)表示累积离差：

式中1≤ t ≤τ 。于是极差 R(τ)定义为：

标准差 S(τ)定义为：

如果采用标准差除极差，相当于将极差“标准化”，消除量纲的影响，于是得到极差与标准差的比值

经过长期的理论分析和模拟试验研究， Hurst 发现如下经验标度关系

式中 H 称为 Hurst 指数。当随机序列{ξ(t) }的前后变化无关、方差有限时， Hurst 等曾经从理论上导出如下关系式：

即有 H=1/2。 Hurst 发现，许多自然现象如江河流量、泥浆沉积等，都有 H>1/2。不过上述规律总是表现在一定的时间尺度范围之内，即存在一个时滞尺度的低限 τmin。基于 Monto Carlo 模拟的结果表明：当 τ>τmin，能够很好地满足 Hurst 经验定律；对于独立的 Gaussian 过程， τmin 变得很大，即需要很大的样本才能进行计算。

2 Hurst指数的Python计算

2.1 code

Hurst指数的Python代码实现如下，本代码是几年前写的，现在我已经看不懂了。

import numpy as np
from osgeo import gdal
import os
"""
计算hurst指数
"""def s(inputdata):# 输入numpy数组n = inputdata.shape[0]t = 0for i in np.arange(n):if i <= (n - 1):for j in np.arange(i + 1, n):if inputdata[j] > inputdata[i]:t = t + 1elif inputdata[j] < inputdata[i]:t = t - 1else:t = treturn tdef beta(inputdata):n = inputdata.shape[0]t = []for i in np.arange(n):if i <= (n - 1):for j in np.arange(i + 1, n):t.append((inputdata[j] - inputdata[i]) / ((j - i) * 1.0))return np.median(t)def Hurst(x):# x为numpy数组n = x.shape[0]t = np.zeros(n - 1)  # t为时间序列的差分for i in range(n - 1):t[i] = x[i + 1] - x[i]mt = np.zeros(n - 1)  # mt为均值序列,i为索引,i+1表示序列从1开始for i in range(n - 1):mt[i] = np.sum(t[0:i + 1]) / (i + 1)# Step3累积离差和极差,r为极差r = []for i in np.arange(1, n):  # i为taocha = []for j in np.arange(1, i + 1):if i == 1:cha.append(t[j - 1] - mt[i - 1])if i > 1:if j == 1:cha.append(t[j - 1] - mt[i - 1])if j > 1:cha.append(cha[j - 2] + t[j - 1] - mt[i - 1])r.append(np.max(cha) - np.min(cha))s = []for i in np.arange(1, n):ss = []for j in np.arange(1, i + 1):ss.append((t[j - 1] - mt[i - 1]) ** 2)s.append(np.sqrt(np.sum(ss) / i))r = np.array(r)s = np.array(s)xdata = np.log(np.arange(2, n))ydata = np.log(r[1:] / s[1:])h, b = np.polyfit(xdata, ydata, 1)return hif __name__ == '__main__':x = np.array([1.59, 1.57, 1.56, 1.54, 1.52, 1.50, 1.47, 1.43, 1.41, 1.40, 1.39])print(Hurst(x))# 0.7486779334192257

2.2 验证

《地理数学方法：基础和应用》一书的【21.2 计算步骤与实例分析】章节中给出了实例分析，使用了中国人均耕地面积变化数据（ 1996－ 2006），将书中给出的数据输入代码中进行计算，得到的结果是0.7486779。书中的计算结果为：0.7487。二者的结果是一致的（保留四位小数），证明了本文代码的正确性。

3 栅格的Hurst指数

有时候可能会需要对时间序列的降水、温度栅格计算Husrt指数图，即对每个像素计算。

以2000-2020年的降水为例，首先将数据合并成一个tif文件，即一个时间对应一个波段，2000年降水量为第1波段，2001年为第2波段，以此类推，2020年为第21波段。

编写一个函数实现影像读写和计算。

def ImageHurst(imgpath,  outtif):"""计算影像的hurst指数:param imgpath: 影像1，多波段:param outtif: 输出结果路径:return: None"""# 读取影像1的信息和数据ds1 = gdal.Open(imgpath)projinfo = ds1.GetProjection()geotransform = ds1.GetGeoTransform()rows = ds1.RasterYSizecolmns = ds1.RasterXSizedata1 = ds1.ReadAsArray()print(data1.shape)src_nodta = ds1.GetRasterBand(1).GetNoDataValue()# 创建输出图像format = "GTiff"driver = gdal.GetDriverByName(format)dst_ds = driver.Create(outtif, colmns, rows, 1,gdal.GDT_Float32)dst_ds.SetGeoTransform(geotransform)dst_ds.SetProjection(projinfo)# 删除对象ds1 = None# 开始计算指数band1 = data1[0]out = band1 * 0 - 2222for row in tqdm(range(rows)):for col in range(colmns):if src_nodta is None:x = data1[:, row, col]hindex  =  Hurst(x)out[row, col] = hindexelse:if band1[row, col] != src_nodta:x = data1[:, row, col]hindex = Hurst(x)out[row, col] = hindex# 写出图像dst_ds.GetRasterBand(1).WriteArray(out)# 设置nodatadst_ds.GetRasterBand(1).SetNoDataValue(-2222)dst_ds = None

4 资料获取

本文完整代码获取方法：公众号【老王GIS】回复关键词【hurst】获取链接。

5 结束语

好朋友们，以上就是今天分享的全部内容了，如有疑问和错误欢迎批评指正。

重标极差分析 Hurst指数计算相关推荐

R语言 Hurst指数计算
GPS filenum=c(11:15,19:24) roads=c("鞍山西道","白堤路","保山道","复康路", ...
Python的Mann-Kendall非参数检验和计算Hurst指数
Mann-Kendall 检验法简称为 M-K 法, 是一种非参数统计检验方法, 可适用于不具有正态分布特征变量的趋势分析[38].假定X1,X2,...Xn为时间序列变量[1],n为时间序列的长度, ...
获取铁矿石和螺纹钢期货数据。对收益率序列进行描述性统计、jb检验，反正是否符合分形市场假说。计算Hurst指数，制定跨品种套利策略，并进行回测，对跨品种套利效果进行评估。寻求改进空间。
源码已上传至github 项目简介获取铁矿石和螺纹钢期货数据.对收益率序列进行描述性统计.jb检验,反正是否符合分形市场假说.计算Hurst指数,制定跨品种套利策略,并进行回测,对跨品种套利效果进行 ...
时间序列中Hurst指数的计算（python代码）
在做时间序列分析时,需要计算Hurst指数,由于Hurst指数计算比较复杂,刚开始懒得自己写,就在github上进行搜索,多是这个代码: from numpy import std, subtract ...
股指的趋势持续研究（Hurst指数）
只贴基本的适合小白的Matlab实现代码,深入的研究除了需要改进算法,我建议好好研究一下混沌与分形,不说让你抓住趋势,至少不会大亏,这个资金盈亏回调我以前研究过. function [line_H,R ...
使用matlab计算hurst指数的代码
您可以使用以下代码来计算Hurst指数: % 加载数据 data = load('your_data.txt');% 计算数据的长度 N = length(data);% 初始化矩阵 rs = zer ...
dfa matlab用法,关于使用MF-DFA方法计算广义Hurst指数的MATLAB操作问题
我在论坛上复制了一个代码,是使用MF-DFA方法计算广义Hurst指数的,但不知道需填入的各个变量的名称,我是零基础,但任务时间很紧,来不及现学,所以想先用来算个数,请各位高手指教,不胜感激! 请问括 ...
泰尔指数计算的stata代码，详细教学，包括相关文献讲解与结果分析
泰尔指数计算的stata代码,详细教学,包括相关文献讲解与结果分析泰尔指数计算详细教学,里面包含文献讲解.结果分析与解读.stata软件操作,使用本人制作的新命令taier11,简单好用,可以计算四种 ...
Hurst指数以及MF-DFA
转:https://uqer.io/home/ https://uqer.io/community/share/564c3bc2f9f06c4446b48393 写在前面 9月的时候说想把arch包加 ...
pythonencoding etf-8_Python 量化分析ETF指数基金投资
Python 量化分析ETF指数基金. 标签(空格分隔): python 量化 ETF tushare pandas 文章目录 Python 量化分析ETF指数基金. 数据获取数据分析在喜马拉雅上 ...

重标极差分析 Hurst指数计算