Eu desenvolvi esses programas há muito tempo atrás em Matlab para calcular o R/S e o V/S da série inteira (Long Range Dependence R/S, V/S, GHE, DFA, DMA). Como muita gente me pede esses códigos estou disponibilizando aqui.
Referências:
Daniel O. Cajueiro, Benjamin M. Tabak. Roberto Andrade Fluctuations in interbank network dynamics. Physical Review. E, Statistical, Nonlinear and Soft Matter Physics, v. 79, p. 037101, 2009.
Daniel O. Cajueiro, Benjamin M. Tabak. The rescaled variance statistic and the determination of the Hurst exponent. Mathematics and Computers in Simulation, Volume 70, Issue 3, 8 November 2005, Pages 172-179
Daniel O. Cajueiro, Benjamin M. Tabak Long-range dependence and market structure. Chaos, Solitons & Fractals, Volume 31, Issue 4, February 2007, Pages 995-1000.
Daniel O. Cajueiro, Benjamin M. Tabak Possible causes of long-range dependence in the Brazilian stock market. Physica A: Statistical Mechanics and its Applications, Volume 345, Issues 3–4, 15 January 2005, Pages 635-645
Daniel O. Cajueiro, Benjamin M. Tabak. Fluctuation dynamics in US interest rates and the role of monetary policy. Finance Research Letters, Volume 7, Issue 3, September 2010, Pages 163-169
Main code
% Programa that evaluates LRD using several different methods (with and without shuffle for R/S and V/S):
% R/S
% V/S
% GHE
% DFA
% DMA
% Input: Series of prices in Excel
% Output: Hurst Exponent
clear all
close all
dados=xlsread('acoes.xls'); % Input your file in Excel
[nLinhas,nColunas]=size(dados);
escalaInicial=5;
escalaFinal=120;
aumentoEscala=1.2;
deltaEscala=5;
numeroMinimoDePontos=10;
blockShuf=10;
for iSerie=1:nColunas
iSerie
% Limpando os dados dos NAN
cont=0;
for i=1:nLinhas
if isnan(dados(i,iSerie))
break;
else
cont=cont+1;
end
end
tamanhoSerie=cont;
serie=log(dados(1:tamanhoSerie,iSerie));
[Hrs,desvio]=frs(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie);
H(1,iSerie)=Hrs;
Hdesvio(1,iSerie)=desvio;
[Hrss,desvio]=frss(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,blockShuf,numeroMinimoDePontos,tamanhoSerie,serie);
H(2,iSerie)=Hrss;
Hdesvio(2,iSerie)=desvio;
[Hvs,desvio]=fvs(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie);
H(3,iSerie)=Hvs;
Hdesvio(3,iSerie)=desvio;
[Hvss,desvio]=fvss(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,blockShuf,numeroMinimoDePontos,tamanhoSerie,serie);
H(4,iSerie)=Hvss;
Hdesvio(4,iSerie)=desvio;
[Hghe,desvio]=fghe(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie);
H(5,iSerie)=Hghe;
Hdesvio(5,iSerie)=desvio;
[Hdfa,desvio]=fdfa(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie);
H(6,iSerie)=Hdfa;
Hdesvio(6,iSerie)=desvio;
[Hdma,desvio]=fdma(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie);
H(7,iSerie)=Hdma;
Hdesvio(7,iSerie)=desvio;
end
save longRangeDependence.dat H -ascII
frs.m:
function [Hrs,desvio]=frs(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie)
for i=2:tamanhoSerie
rserie(i,1)=serie(i,1)-serie(i-1,1);
end
tamanhoSerie=tamanhoSerie-1;
clear serie;
serie=rserie;
clear rserie;
% Zerando os vetores
sigmaRS=[];
escalas=[];
% Calculo do DFA
escala=escalaInicial;
numeroEscalas=0;
while ((escala<tamanhoSerie-numeroMinimoDePontos)&(escala<=escalaFinal))
numeroEscalas=numeroEscalas+1;
sigmaRS(numeroEscalas,1)=0;
for i=1:tamanhoSerie-escala-1
%Definicao da serie dentro da escala
serieEscala=serie(i:i+escala-1,1);
m=mean(serieEscala);
s=sqrt(var(serieEscala));
x=[];
x(1,1)=serieEscala(1,1)-m; % (j*tau) indica em que escala esta se iniciando a contagem
% De fato, acima se esta calculando o R
for j=2:escala
x(j,1)=x(j-1,1)+serieEscala(j)-m;
end
% Calculo de max e minx dentro da escala
maxx=max(x);
minx=min(x);
r=maxx-minx;
sigmaRS(numeroEscalas,1)=sigmaRS(numeroEscalas,1)+(r/s);
end
sigmaRS(numeroEscalas,1)=sigmaRS(numeroEscalas,1)/(tamanhoSerie-escala+1);
escalas(numeroEscalas,1)=escala;
if (escala*aumentoEscala-escala>deltaEscala)
escala=round(escala*aumentoEscala);
else
escala=escala+deltaEscala;
end
end
X=[ones(numeroEscalas,1) log(escalas)];
Y=log(sigmaRS);
[B,Bint] = regress(Y,X);
Hrs=B(2);
desvio=Bint(2,2)-B(2);
frss.m
function [Hrss,desvio]=frss(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,block_shuf,numeroMinimoDePontos,tamanhoSerie,serie)
for i=2:tamanhoSerie
rserie(i-1,1)=serie(i,1)-serie(i-1,1);
end
tamanhoSerie=tamanhoSerie-1;
clear serie;
serie=rserie;
clear rserie;
% Shuffling data
conts=0;
for i=1:block_shuf:tamanhoSerie-block_shuf
aux=[];
conts=0;
for j=i:i+block_shuf-1
conts=conts+1;
aux(conts)=serie(j,1);
end
auxrand=randperm(block_shuf);
conts=0;
for j=i:i+block_shuf-1
conts=conts+1;
serie(j,1)=aux(auxrand(conts));
end
end
% Zerando os vetores
sigmaRS=[];
escalas=[];
% Calculo do DFA
escala=escalaInicial;
numeroEscalas=0;
while ((escala<tamanhoSerie-numeroMinimoDePontos)&(escala<=escalaFinal))
numeroEscalas=numeroEscalas+1;
sigmaRS(numeroEscalas,1)=0;
for i=1:tamanhoSerie-escala-1
%Definicao da serie dentro da escala
serieEscala=serie(i:i+escala-1,1);
m=mean(serieEscala);
s=sqrt(var(serieEscala));
x=[];
x(1,1)=serieEscala(1,1)-m; % (j*tau) indica em que escala esta se iniciando a contagem
% De fato, acima se esta calculando o R
for j=2:escala
x(j,1)=x(j-1,1)+serieEscala(j)-m;
end
% Calculo de max e minx dentro da escala
maxx=max(x);
minx=min(x);
r=maxx-minx;
sigmaRS(numeroEscalas,1)=sigmaRS(numeroEscalas,1)+(r/s);
end
sigmaRS(numeroEscalas,1)=sigmaRS(numeroEscalas,1)/(tamanhoSerie-escala+1);
escalas(numeroEscalas,1)=escala;
if (escala*aumentoEscala-escala>deltaEscala)
escala=round(escala*aumentoEscala);
else
escala=escala+deltaEscala;
end
end
X=[ones(numeroEscalas,1) log(escalas)];
Y=log(sigmaRS);
[B,Bint] = regress(Y,X);
Hrss=B(2);
desvio=Bint(2,2)-B(2);
fvs.m
function [Hvs,desvio]=fvs(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie)
for i=2:tamanhoSerie
rserie(i-1,1)=serie(i,1)-serie(i-1,1);
end
tamanhoSerie=tamanhoSerie-1;
clear serie;
serie=rserie;
clear rserie;
% Zerando os vetores
sigma2VS=[];
escalas=[];
% Calculo do VS
escala=escalaInicial;
numeroEscalas=0;
while ((escala<tamanhoSerie-numeroMinimoDePontos)&(escala<=escalaFinal))
numeroEscalas=numeroEscalas+1;
sigma2VS(numeroEscalas,1)=0;
for i=1:tamanhoSerie-escala-1
%Definicao da serie dentro da escala
serieEscala=serie(i:i+escala-1,1);
m=mean(serieEscala);
s=sqrt(var(serieEscala));
sumk1=0;
sumk2=0;
for j=1:escala
suml=0;
for k=1:j
suml=suml+serieEscala(k)-m;
end
sumk1=sumk1+suml^2;
sumk2=sumk2+suml;
end
sigma2VS(numeroEscalas)=sigma2VS(numeroEscalas)+(1/((s^2)*escala))*(sumk1-(1/escala)*((sumk2)^2));
end
sigma2VS(numeroEscalas,1)=sigma2VS(numeroEscalas,1)/(tamanhoSerie-escala-1);
escalas(numeroEscalas,1)=escala;
if (escala*aumentoEscala-escala>deltaEscala)
escala=round(escala*aumentoEscala);
else
escala=escala+deltaEscala;
end
end
X=[ones(numeroEscalas,1) log(escalas)];
Y=0.5*log(sigma2VS);
[B,Bint] = regress(Y,X);
Hvs=B(2);
desvio=Bint(2,2)-B(2);
fvss.m
function [Hvss,desvio]=fvss(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,block_shuf,numeroMinimoDePontos,tamanhoSerie,serie)
for i=2:tamanhoSerie
rserie(i-1,1)=serie(i,1)-serie(i-1,1);
end
tamanhoSerie=tamanhoSerie-1;
clear serie;
serie=rserie;
clear rserie;
% Shuffling data
conts=0;
for i=1:block_shuf:tamanhoSerie-block_shuf
aux=[];
conts=0;
for j=i:i+block_shuf-1
conts=conts+1;
aux(conts)=serie(j,1);
end
auxrand=randperm(block_shuf);
conts=0;
for j=i:i+block_shuf-1
conts=conts+1;
serie(j,1)=aux(auxrand(conts));
end
end
% Zerando os vetores
sigma2VS=[];
escalas=[];
% Calculo do VS
escala=escalaInicial;
numeroEscalas=0;
while ((escala<tamanhoSerie-numeroMinimoDePontos)&(escala<=escalaFinal))
numeroEscalas=numeroEscalas+1;
sigma2VS(numeroEscalas,1)=0;
for i=1:tamanhoSerie-escala-1
%Definicao da serie dentro da escala
serieEscala=serie(i:i+escala-1,1);
m=mean(serieEscala);
s=sqrt(var(serieEscala));
sumk1=0;
sumk2=0;
for j=1:escala
suml=0;
for k=1:j
suml=suml+serieEscala(k)-m;
end
sumk1=sumk1+suml^2;
sumk2=sumk2+suml;
end
sigma2VS(numeroEscalas)=sigma2VS(numeroEscalas)+(1/((s^2)*escala))*(sumk1-(1/escala)*((sumk2)^2));
end
sigma2VS(numeroEscalas,1)=sigma2VS(numeroEscalas,1)/(tamanhoSerie-escala-1);
escalas(numeroEscalas,1)=escala;
if (escala*aumentoEscala-escala>deltaEscala)
escala=round(escala*aumentoEscala);
else
escala=escala+deltaEscala;
end
end
X=[ones(numeroEscalas,1) log(escalas)];
Y=0.5*log(sigma2VS);
[B,Bint] = regress(Y,X);
Hvss=B(2);
desvio=Bint(2,2)-B(2);
fghe.m
function [Hghe,desvio]=fghe(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie)
% Zerando os vetores
K2=[];
escalas=[];
% Calculo do DFA
escala=escalaInicial;
numeroEscalas=0;
while ((escala<tamanhoSerie-numeroMinimoDePontos)&(escala<=escalaFinal))
numeroEscalas=numeroEscalas+1;
sigma2fa(numeroEscalas,1)=0;
soma1=0;
soma2=0;
for i=1:tamanhoSerie-escala-1
soma1=soma1+(abs(serie(i+escala)-serie(i)))^2;
soma2=soma2+(abs(serie(i)))^2;
end
soma1=soma1/(tamanhoSerie-escala-1);
soma2=soma2/(tamanhoSerie-escala-1);
K2(numeroEscalas,1)=soma1/soma2;
escalas(numeroEscalas,1)=escala;
if (escala*aumentoEscala-escala>deltaEscala)
escala=round(escala*aumentoEscala);
else
escala=escala+deltaEscala;
end
end
X=[ones(numeroEscalas,1) log(escalas)];
Y=0.5*log(K2);
[B,Bint] = regress(Y,X);
Hghe=B(2);
desvio=Bint(2,2)-B(2);
fdfa.m:
function [Hdfa,desvio]=fdfa(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie)
seriex=(1:tamanhoSerie)';
% Zerando os vetores
sigma2fa=[];
escalas=[];
% Calculo do DFA
escala=escalaInicial;
numeroEscalas=0;
while ((escala<tamanhoSerie-numeroMinimoDePontos)&(escala<=escalaFinal))
numeroEscalas=numeroEscalas+1;
sigma2fa(numeroEscalas,1)=0;
for i=1:tamanhoSerie-escala-1
%Definicao da serie dentro da escala
serieEscala=serie(i:i+escala-1,1);
serieEscalax=seriex(i:i+escala-1,1);
% Calculo da regressao dentro da escala
regressX=[ones(escala,1) serieEscalax];
regressY=serieEscala;
coeficiente=regress(regressY,regressX);
% Calculo do dfa dentro de cada escala
for j=1:escala
sigma2fa(numeroEscalas,1)=sigma2fa(numeroEscalas,1)+((serie(i+j-1)-coeficiente(1)-coeficiente(2)*seriex(i+j-1))^2)/escala; %trocar abs
end
end
sigma2fa(numeroEscalas,1)=sigma2fa(numeroEscalas,1)/(tamanhoSerie-escala-1);
escalas(numeroEscalas,1)=escala;
if (escala*aumentoEscala-escala>deltaEscala)
escala=round(escala*aumentoEscala);
else
escala=escala+deltaEscala;
end
end
X=[ones(numeroEscalas,1) log(escalas)];
Y=0.5*log(sigma2fa);
[B,Bint] = regress(Y,X);
Hdfa=B(2);
desvio=Bint(2,2)-B(2);
fdma:
function [Hdma,desvio]=fdma(escalaInicial,escalaFinal,aumentoEscala,deltaEscala,numeroMinimoDePontos,tamanhoSerie,serie)
% Zerando os vetores
sigma2MA=[];
escalas=[];
% Calculo do DMA
escala=escalaInicial;
numeroEscalas=0;
while ((escala<tamanhoSerie-numeroMinimoDePontos)&(escala<=escalaFinal))
numeroEscalas=numeroEscalas+1;
sigma2MA(numeroEscalas,1)=0;
for i=escala+1:tamanhoSerie
%Definicao da serie dentro da escala
serieEscala=serie(i-escala:i-1,1);
% Calculo da media movel dentro da escala
movingAverage=0;
for j=1:escala
movingAverage=movingAverage+serieEscala(j);
end
movingAverage=movingAverage/escala;
% Calculo do Sigma2MA dentro de cada escala
sigma2MA(numeroEscalas,1)=sigma2MA(numeroEscalas,1)+(serie(i)-movingAverage)^2;
end
sigma2MA(numeroEscalas,1)=sigma2MA(numeroEscalas,1)/(tamanhoSerie-escala-1);
escalas(numeroEscalas,1)=escala;
if (escala*aumentoEscala-escala>deltaEscala)
escala=round(escala*aumentoEscala);
else
escala=escala+deltaEscala;
end
end
X=[ones(numeroEscalas,1) log(escalas)];
Y=0.5*log(sigma2MA);
[B,Bint] = regress(Y,X);
Hdma=B(2);
desvio=Bint(2,2)-B(2);