Hurst指数估计方法（时域）—

求Hurst指数的意义：
作为判断时间序列数据遵从随机游走还是有偏的随机游走过程的指标。简单来说就是为了通过计算的这个指数判断一下未来的趋势，也有用来判断时间序列平稳性的。
Hurst指数估计的方法有很多，按照时域和频域可以分为如下几个【1】：
时域：方差-时间法、聚合序列绝对值法、R/S法（用得最多）和趋势波动分析（DFA）
频域：周期图法、Whittle估计法和小波分析法
Hurst指数在不同范围下的意义：【2】

DFA的计算原理和计算步骤（关于DFA的计算版本有很多，这里只是其中的一种）
已知：一个时间序列，长度为N，记作X，时间序列表示为{xi}，i = 1，2，…，N；
①取序列，选择合适的区间长度，也就是时间窗s，将序列划分为长度为s的不重叠等长度子区间，长度为N的子序列就会被分成Ns = N/s个段，但是想一下极端的事件，时间序列的个数是个质数怎么办，不要后面除不尽的又好浪费数据，于是为保证序列信息不丢失，可以取两次数据，第一次丢掉最后几个取不到的数据，第二次就倒着看，正着算，丢掉前面几个取不到的数据。
Ns=floor(N/s),其中floor()表示向下取整Ns = floor(N/s),其中floor()表示向下取整 Ns=floor(N/s),其中floor()表示向下取整
m1=Ns∗sm1 = Ns*s m1=Ns∗s
X1=x1,x2,...xm1X1 = {x1,x2,...xm1} X1=x1,x2,...xm1
m2=N−m1m2 = N-m1 m2=N−m1
X2=xm2,....xNX2 = {xm2,....xN} X2=xm2,....xN
②分别对新序列X1和X2计算累积离差，生成新的序列Y1和Y2

③对每个子区间进行多项式拟合，得到局部趋势韩式，并计算消除趋势的序列：

Python代码

# author: Dominik Krzeminski (dokato)import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as ss
import xlrd# detrended fluctuation analysisdef calc_rms(x, scale):"""windowed Root Mean Square (RMS) with linear detrending.Args:-----*x* : numpy.arrayone dimensional data vector*scale* : intlength of the window in which RMS will be calculaedReturns:--------*rms* : numpy.arrayRMS data in each window with length len(x)//scale"""# making an array with data divided in windowsshape = (x.shape[0]//scale, scale)X = np.lib.stride_tricks.as_strided(x,shape=shape)# vector of x-axis points to regressionscale_ax = np.arange(scale)rms = np.zeros(X.shape[0])for e, xcut in enumerate(X):coeff = np.polyfit(scale_ax, xcut, 1)xfit = np.polyval(coeff, scale_ax)# detrending and computing RMS of each windowrms[e] = np.sqrt(np.mean((xcut-xfit)**2))return rmsdef dfa(x, scale_lim=[4,10], scale_dens=0.25, show=False):"""Detrended Fluctuation Analysis - measures power law scaling coefficientof the given signal *x*.More details about the algorithm you can find e.g. here:Hardstone, R. et al. Detrended fluctuation analysis: A scale-free view on neuronal oscillations, (2012).Args:-----*x* : numpy.arrayone dimensional data vector*scale_lim* = [5,9] : list of length 2 boundaries of the scale, where scale means windows among which RMSis calculated. Numbers from list are exponents of 2 to the powerof X, eg. [5,9] is in fact [2**5, 2**9].You can think of it that if your signal is sampled with F_s = 128 Hz,then the lowest considered scale would be 2**5/128 = 32/128 = 0.25,so 250 ms.*scale_dens* = 0.25 : floatdensity of scale divisions, eg. for 0.25 we get 2**[5, 5.25, 5.5, ... ] *show* = Falseif True it shows matplotlib log-log plot.Returns:--------*scales* : numpy.arrayvector of scales (x axis)*fluct* : numpy.arrayfluctuation function values (y axis)*alpha* : floatestimation of DFA exponent"""# cumulative sum of data with substracted offsety = np.cumsum(x - np.mean(x))scales = (2**np.arange(scale_lim[0], scale_lim[1], scale_dens)).astype(np.int)# scales = list(range(4,17))fluct = np.zeros(len(scales))# computing RMS for each windowfor e, sc in enumerate(scales):fluct[e] = np.sqrt(np.mean(calc_rms(y, sc)**2))# fitting a line to rms datacoeff = np.polyfit(np.log2(scales), np.log2(fluct), 1)if show:fluctfit = 2**np.polyval(coeff,np.log2(scales))plt.loglog(scales, fluct, 'bo')plt.loglog(scales, fluctfit, 'r', label=r'$\alpha$ = %0.2f'%coeff[0])plt.title('DFA')plt.xlabel(r'$\log_{10}$(time window)')plt.ylabel(r'$\log_{10}$<F(t)>')plt.legend()plt.show()return scales, fluct, coeff[0]if __name__=='__main__':# n = 1000# x = np.random.randn(n)# # computing DFA of signal envelope# x = np.abs(ss.hilbert(x))# 径流输入YC_flow = 'SX-data/YC.xlsx'input_filePath = YC_flowinput_data = []wb = xlrd.open_workbook(filename=input_filePath)print(wb.sheet_names())for year in wb.sheet_names():# year_data = []sheet = wb.sheet_by_name(year)cols = sheet.col_values(1)for col in cols:# year_data.append(col)input_data.append(col)scales, fluct, alpha = dfa(input_data, show=1)print(scales)print(fluct)print("DFA exponent: {}".format(alpha))

matlab代码

function[A,F] = DFA_fun(data,pts,order)% -----------------------------------------------------
% DESCRIPTION:
% Function for the DFA analysis.% INPUTS:
% data: a one-dimensional data vector.
% pts: sizes of the windows/bins at which to evaluate the fluctuation
% order: (optional) order of the polynomial for the local trend correction.
% if not specified, order == 1;% OUTPUTS:
% A: a 2x1 vector. A(1) is the scaling coefficient "alpha",
% A(2) the intercept of the log-log regression, useful for plotting (see examples).
% F: A vector of size Nx1 containing the fluctuations corresponding to the
% windows specified in entries in pts.
% -----------------------------------------------------% Checking the inputs
if nargin == 2order = 1;
endsz = size(data);
if sz(1)< sz(2)data = data';
endexit = 0;if min(pts) == order+1disp(['WARNING: The smallest window size is ' num2str(min(pts)) '. DFA order is ' num2str(order) '.'])disp('This severly affects the estimate of the scaling coefficient')disp('(If order == [] (so 1), the corresponding fluctuation is zero.)')
elseif min(pts) < (order+1)disp(['ERROR: The smallest window size is ' num2str(min(pts)) '. DFA order is ' num2str(order) ':'])disp(['Aborting. The smallest window size should be of ' num2str(order+1) ' points at least.'])exit = 1;
endif exit == 1return
end% DFA
npts = numel(pts); %返回数组pts中元素个数F = zeros(npts,1);
N = length(data);for h = 1:nptsw = pts(h);n = floor(N/w);Nfloor = n*pts(h);D = data(1:Nfloor);y = cumsum(D-mean(D));bin = 0:w:(Nfloor-1);vec = 1:w;coeff = arrayfun(@(j) polyfit(vec',y(bin(j) + vec),order),1:n,'uni',0);y_hat = cell2mat(cellfun(@(y) polyval(y,vec),coeff,'uni',0));F(h)  = mean((y - y_hat').^2)^0.5;endA = polyfit(log(pts),log(F)',1);end

[1]刘付斌, 高相铭. 基于EEMD与DFA的Hurst指数估计[J]. 测控技术, 2013(10):101-104.
[2]袁全勇,杨阳,李春,阚威,叶柯华.基于Hurst指数的风速时间序列研究[J].应用数学和力学,2018,39(07):798-810.

Hurst指数估计方法（时域）——DFA相关推荐

dfa matlab用法,关于使用MF-DFA方法计算广义Hurst指数的MATLAB操作问题
我在论坛上复制了一个代码,是使用MF-DFA方法计算广义Hurst指数的,但不知道需填入的各个变量的名称,我是零基础,但任务时间很紧,来不及现学,所以想先用来算个数,请各位高手指教,不胜感激! 请问括 ...
matlab 双边沿滤波,一种基于数字PWM发生器的左增长双边沿UPWM信号频谱估计方法与流程...
本发明涉及数字D类音频功放领域,尤其涉及一种由数字音频信号调制得到的左增长双边沿均匀采样脉冲宽度调制信号的频谱估计方法. 背景技术: 数字D类音频功放的电源效率相比A类.B类和AB类等线性音频功放较高 ...
Hurst指数以及MF-DFA
转:https://uqer.io/home/ https://uqer.io/community/share/564c3bc2f9f06c4446b48393 写在前面 9月的时候说想把arch包加 ...
Python的Mann-Kendall非参数检验和计算Hurst指数
Mann-Kendall 检验法简称为 M-K 法, 是一种非参数统计检验方法, 可适用于不具有正态分布特征变量的趋势分析[38].假定X1,X2,...Xn为时间序列变量[1],n为时间序列的长度, ...
指数平滑方法(一次指数平滑、二次指数平滑、三次指数平滑)：理论、代码、参数介绍（全）
@创建于:20210324 @修改于:20210324 文章目录特别说明参考来源包版本号 1.简介 2.一次指数平滑 2.1 理论介绍 2.2 代码展示 2.3 参数介绍 3. 二次指数平滑 3 ...
指数平滑方法深度解析（一次二次三次）
简书同步参考链接指数平滑方法说起来感觉挺简单的,不就是几期求均值吗,但是你知道在Eviews里做指数平滑模型的时候,1.他的初始值是如何确定的吗?2.初始值的确定方法可以按照我们想的去改变吗? 3. ...
两发两收信道估计方法
本文用于记录2发2收系统的信道估计方法,采用最基础的LS估计算法,主要是记录多发多收系统该如何处理多信道的估计汇总. 如图所示,是一个2发2收系统的基本模型从频率域分析:其中发射端Tx1T_{x1} ...
信号频谱质心matlab,对称频谱信号的中心频率的质心估计方法
对称频谱信号的中心频率的质心估计方法 [专利摘要]对称频谱信号的中心频率的质心估计方法.由于能量泄漏和栅栏效应的影响使得用FFT得到的离散频谱直接估计信号的中心频率会产生较大误差.本发明方法包括:步骤 ...
广义Hurst指数与分维D关系的范例解释
金融领域常见的Hurst指数实指广义Hurst指数(Generalized Hurst Exponent),由Mandelbrot提出.与Hurst本人最早提出的Hurst指数并不相同. 一.Hurs ...
基于训练符号的频偏估计方法（SC-FDE/OFDM）
1.设置归一化频偏时域(使用l两段完全一样的chu序列,这里设置每一段长度为1024),基带传输 x(k)=x(k+N),k=1,2,...,1024x(k )=x(k+N),k=1,2,...,1 ...

Hurst指数估计方法（时域）——DFA

Hurst指数估计方法（时域）——DFA相关推荐

最新文章

热门文章