200字范文 > matlab DFA算法计算Hurst指数

matlab DFA算法计算Hurst指数

时间：2018-09-23 11:00:16

Hurst 指数是分形市场理论中最关键的指标，被用于描述时间序列的长记忆性程度，当H值等于0.5则价格序列表现为随机游走，当大于0.5则表明时间序列具有长记忆性，小于0.5则表明存在均值回复特性。

算法描述

计算 Hurst 指数有多种方法，本文采用 matlab 实现DFA(Detrended FluctuationAnalysis) 即降趋脉动分析法。该方法描述如下：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CcvZyh4r-1586690321723)(http://p0yg4u41b./-10-30-22-33-52.png)]

运用该方法求解以下问题：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-YVWBisdQ-1586690321724)(http://p0yg4u41b./-10-30-22-36-06.png)]

算法实现

核心算法

通过分析算法，我们可以知道需要用s长度的窗口去切分数据。因此我们先将窗口s作为一个固定的输入参数，求解出其对应的 F(s) 。编写 DFA_core函数，该函数代码如下;

function f = DFA_core(data,s,order)% data 数据，s 窗口长度，order 所估计的多项式次数%% 构造新时间序列temp = data - mean(data);data1 = cumsum(temp);% 原时间序列长len = length(data);% 共可分为n组n = floor(len/s);% 真正用到的时间序列长len1 = n * s;%% 函数拟合% 预分配拟合参数fitcoef = nan(n,2);% 拟合x = (1:s)';for i = 1:ny = data1((i-1)*s+1:i*s);fitcoef(i,:) = polyfit(x,y,order);end%% 计算每个子序列降趋后剩余序列的方差(按题中的公式表述为方差应该不正确）fk = nan(n,1);for i = 1:n% 计算拟合多项式的值g = polyval(fitcoef(i,:),x);y = data1((i-1)*s+1:i*s);temp = y-g;fk(i) = mean(temp .* temp);end%% 计算样本的总体降趋脉动的均方根f2 = mean(fk);f= sqrt(f2);end

Hurst 指数估计

通过调用上文中编写的 DFA_core 函数，可以得到 s 序列对应的 F(s) 序列的值。我们知道他们是呈现出幂函数关系，所以对其取对数，有logF(s) = Hlog(s) + log(a)。代码如下：

function H = myDFA(data,s,order)% order 多项式阶数n = length(s);f = nan(n,1);% 调用DFA_core计算不同s对应的脉动ffor i = 1:nf(i) = DFA_core(data,s(i),order);endtoc%% 拟合出H指数% 将幂函数转换成对数函数logf = log(f);logs = log(s');% h = polyfit(logs,logf,1);% H = h(1);X=[ones(n,1) logs];y=logf;b=regress(y,X);H = b(2);%% plotfigure;plot(logs,logf,'*');xlabel('logs');ylabel('logf');hold onY = H*logs + b(1);plot(logs,Y);legend("散点","拟合曲线");end

问题求解

编写脚本，解决问题：

%% 问题二脚本% load dataclose allclearclcload stock_dataticClose = data(:,end);% 多项式阶数为1order = 1;% s为窗口长度，题目中没有给定范围，注意窗口必须大于多% 项式阶数，否则无法估计。暂选 10：1000s = 10:1500;%% % 计算日度对数收益率r = price2ret(Close);H_r = myDFA(r,s,order);%% 计算日度波动率Open = data(:,2);High = data(:,3);Low = data(:,4);ht = High - Open;lt = Low - Open;ct = Close - Open;k1 = 0.511;k2 = 0.019;k3 = 0.383;vt = k1 * ((ht-lt) .^2) -k2 * (ct .* (ht + lt) - 2*ht .* lt) - k3 * (ct .^2);H_v = myDFA(vt,s,order);%% 输出fprintf('The Hurst index of returns is %.3f\n',H_r)fprintf('The Hurst index of volatilities is %.3f\n',H_v)toc