MultiFractal Detrended Fluctuation Analysis: parallel version

MFDFA2.R

library(MFDFA)
library(foreach)
library(doParallel)

MFDFA2<-function (tsx, scale, m = 1, q) {
  
  poly_fit.val<-function (x, y, n) {
    formule <- lm(as.formula(paste("y~", paste("I(x^", 1:n, ")", 
                                               sep = "", collapse = "+"))))
    res1 <- coef(formule)
    poly.res <- res1[length(res1):1]
    suppressWarnings(res2 <- predict(formule, interval = "prediction"))
    poly.eva <- res2[, 1]
    allres <- list(polyfit = round(poly.res, 4), polyval = poly.eva)
    return(allres)
  }
  
  X <- cumsum(tsx - mean(tsx))
  seg <- c()
  qRMS <- list()
  Fq <- c()
  Fqi <- list()
  Hq <- c()
  qRegLine <- list()
  RMSvi <- list()
  
  cores=detectCores()
  cl <- makeCluster(cores[1]-1) 
  registerDoParallel(cl)
  
  Fqi<-foreach(i=1:length(scale)) %dopar% {
    seg[i] <- floor(length(X)/scale[i])
    rmvi <- c()
    for (vi in 1:seg[i]) {
      Index = ((((vi - 1) * scale[i]) + 1):(vi * scale[i]))
      polyft <- poly_fit.val(Index, X[Index], m)
      C <- polyft$polyfit
      fit <- polyft$polyval
      rmvi[vi] <- sqrt(mean((fit - X[Index])^2))
      RMSvi[[i]] <- rmvi
    }
    for (nq in 1:length(q)) {
      rod <- RMSvi[[i]]^q[nq]
      qRMS[[nq]] <- rod
      Fq[nq] <- mean(qRMS[[nq]])^(1/q[nq])
    }
    if (any(q == 0)) {
      Fq[which(q == 0)] <- exp(0.5 * mean(log(RMSvi[[i]]^2)))
    }
    Fqi[[i]] <- Fq
  }
  
  stopCluster(cl)
  
  
  Fqi <- Reduce("rbind", Fqi)
  for (nq in 1:length(q)) {
    polyft <- poly_fit.val(log2(scale), log2(Fqi[, nq]), 1)
    C <- polyft$polyfit
    Hq[nq] <- C[1]
    qRegLine[[nq]] <- polyft$polyval
  }
  qRegLine <- as.data.frame(Reduce("cbind", qRegLine))
  tq <- Hq * q - 1
  hq <- diff(tq)/(q[2] - q[1])
  Dq <- (q[1:(length(q) - 1)] * hq) - tq[1:(length(tq) - 1)]
  return(list(Hq = Hq, tau_q = tq, spec = data.frame(hq = hq,Dq = Dq), Fqi = Fqi, line = qRegLine))
}

https://github.com/mlaib/MFDFA#little-comparison