MF-DFA from Matalb to Julia implementation

I'm working with the MF-DFA method in Matlab but I need to implement it in Julia. The goal is tho obtain the Hurst exponent of the S&P 500. The matlab code is as follows:

sp500 = readtable('sp500_Nasdaq.csv','PreserveVariableNames', true) ;
spClose = table2array(sp500(:,2))  ;

SP1=cumsum(spClose - mean(spClose)) ;
SP1_ordinary=sqrt(mean(SP1.^2));
X=cumsum(spClose-mean(spClose));
X=transpose(X);
scale=[16,32,64,128,256,512,1024];
q=[-5,-3,-1,0,1,3,5];
m=1;
for ns=1:length(scale),
segments(ns)=floor(length(X)/scale(ns));
for v=1:segments(ns)
  Index=( ( ( (v-1)*scale(ns) )+1):(v*scale(ns)));
  C = polyfit(Index,X(Index),m) ;
  fit=polyval(C,Index);
  RMS{ns}(v)=sqrt(mean((X(Index)-fit).^2));
end
for nq=1:length(q),
  qRMS{nq,ns}=RMS{ns}.^q(nq);
  Fq(nq,ns)=mean(qRMS{nq,ns}).^(1/q(nq));
end
  Fq(q==0,ns)=exp(0.5*mean(log(RMS{ns}.^2)));
end

The Julia the code looks like this:

using DelimitedFiles, TimeSeries, Plots, DelimitedFiles, Plots, StatsBase
using Polynomials, LinearAlgebra, CSV, DataFrames

sp500 = CSV.read("sp500_Nasdaq.csv", DataFrame) 
sp500_V = values(sp500[:,2]) 
SP1 = cumsum(sp500_V .- mean(sp500_V) ) ; 
SP1_Ord = sqrt(mean(SP1.^2)) ;
X = SP1 ;
X = X';
function polyfit(xVals,yVals)
   n = length(xVals)
   xBar, yBar = @fastmath mean(xVals), mean(yVals)
   sXX, sXY = @fastmath ones(n)'*(xVals.-xBar).^2 , dot(xVals.-xBar,yVals.-yBar)
   b1A = @fastmath sXY/sXX
   b0A = @fastmath yBar - b1A*xBar
return b0A, b1A
end

scales = [16,32,64,128,256,512,1024];
q = [-5,-3,-1,0,1,3,5] ;
segments = zeros(Int64, (1,length(scales)))
global qRMS = zeros( length(q) ,length(scales)  ) ;
global Fq = zeros( length(q) , length(scales) ) ;
@inbounds for ns = 1:length(scales)
global segments[ns] = Int(floor( length(X)/scales[ns] ) ) ;  
global Index = Array{UnitRange{Int128}}(undef, (segments[ns], length(scales))  ) ;
global ft = zeros(Float64, (segments[ns], length(scales) ) ) ;
global RMS = zeros(Float64, (length(scales) ,segments[ns] ) ) ;

@inbounds  for v=1:segments[ns]
    global RMSk = Array{Float64}[] ;
    Index =  ( ( (v-1)*scales[ns] ) + 1 ):( v*scales[ns] ) ;
    global C = polyfit( Index, X[Index]) ;
    global p = Polynomial(C)
    ft =p.(Index);
    RMS[ns,v] = sqrt(mean((X[Index] .- ft).^2));
    push!(RMSk,RMS )
end
@inbounds for nq = 1:length(q)
    qRMS[nq,ns] = RMS[ns].^q[nq];
    Fq[nq,ns] = mean( qRMS[nq,ns] ).^(1/q[nq] );
end
Fq[findall(x->x==0, q)[1], ns] = exp( 0.5*mean(log.(RMS[ns].^2) ) ) ;
end

The thing is that the array RMS in Matlab's code is an array of arrays like this:

RMS =

1×7 cell array

{1×159 double}    {1×79 double}    {1×39 double}    {1×19 double}    {1×9 double}    {1×4 double}    {1×2 double}

But Julia returns only the last array

RMS
7×2 Matrix{Float64}:
 0.0      0.0
 0.0      0.0
 0.0      0.0
 0.0      0.0
 0.0      0.0
 0.0      0.0
62178.0  18238.2

How can I obtain the same output as in Matlab? How can you store arrays into arrays in Julia?

1 answer

  • answered 2021-11-23 23:36 utello10

    The solutions for this is to use RMScell = Array{Float64}[] is equivalent to a Matlab cell array

    using DelimitedFiles, TimeSeries, Plots, DelimitedFiles, Plots
    using Polynomials, LinearAlgebra, CSV, DataFrames, StatsBase
    
    sp500 = CSV.read("sp500_Nasdaq.csv", DataFrame) ;
    sp500_V = values(sp500[:,2]) ;
    SP1 = cumsum( sp500_V .- mean(sp500_V) ) ; 
    SP1_Ord = sqrt( mean(SP1.^2) ) ;
    
    X = SP1 ;
    X = X' ;
    function polyfit(xVals,yVals)
       n = length(xVals)
       xBar, yBar = @fastmath mean(xVals), mean(yVals)
       sXX, sXY = @fastmath ones(n)'*(xVals.-xBar).^2 , dot(xVals.-  xBar,yVals.-yBar)
       b1A = @fastmath sXY/sXX
       b0A = @fastmath yBar - b1A*xBar
    return b0A, b1A
    end
    """ Multifractal detrended fluctuation analysis of time series """
    scales = [16,32,64,128,256,512,1024];
    q = [-5,-3,-1,0,1,3,5] ;
    segments = zeros(Int64, (1,length(scales))) ;
    global qRMS = zeros( length(q) ,length(scales)  ) ;
    global Fq = zeros( length(q) , length(scales) ) ;
    global RMScell = Array{Float64}[] ;
    global qRMScell =[] ;
    global segmentsFq = [] ;
    @inbounds for ns = 1:length(scales)
    global segments[ns] = Int(floor( length(X)/scales[ns] ) ) ;  
    global ft = zeros(Float64, (segments[ns], length(scales) ) ) ;
    global RMS = zeros(Float64, segments[ns]);
    @inbounds  for v=1:segments[ns]
        global Index = ( (v-1)*scales[ns] ) + 1: v*scales[ns] ;
        global C = polyfit( Index, X[Index]) ;
        global p = Polynomial(C) ;
        ft =p.(Index) ;
        RMS[v] = sqrt(mean((X[Index] .- ft).^2))  ;
        end
        l = deepcopy(RMS)
        push!(RMScell,l)
        global IndexFq = ((ns-1)*length(q) ) + 1 : ns*length(q) ;
        push!(segmentsFq, IndexFq) ;
    @inbounds for nq = 1:length(q)
        l = RMScell[ns].^q[nq]
        r = deepcopy(l) ;
        push!(qRMScell, r) ;
        end
    @inbounds for nq = 1: length(scales)
        Fq[nq,ns] = mean( qRMScell[segmentsFq[ns]][nq] ).^(1/q[nq] ) ;
        end
    Fq[findall(x->x==0, q)[1], ns] = exp( 0.5*mean(log.(RMScell[ns].^2) ) ) ;
    

    end

    Hq = zeros( Float64,length(q) ) ;
    global qRegLine = Array{Float64}[] ;
    for  nq = 1:length(q)
    global C = polyfit( log2.(scales),log2.(Fq[nq,:]) ) ;
    Hq[nq] = C[2] ;
    global p = Polynomial(C) ;
    push!( qRegLine, p.( log2.(scales) ) )
    end
    
    
    tq = Hq.*q .- 1 ;
    hq = diff(tq)./(q[2]-q[1]) ;
    Dq = ( q[1:end-1].*hq ) - tq[1:end-1] ;
    

How many English words
do you know?
Test your English vocabulary size, and measure
how many words do you know
Online Test
Powered by Examplum