MFDFA_plot.R

##################################
## Suggestion of output plot: ####
##################################

## Supplementary functions: #####
reset <- function(){
par(mfrow=c(1, 1), oma=rep(0, 4), mar=rep(0, 4), new=TRUE)
plot(0:1, 0:1, type="n", xlab="", ylab="", axes=FALSE)}

poly_fit<-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]
  allres<-list(polyfit=poly.res, model1=formule)
  return(allres)}
##################################

#####################
## Output plots: ####
dev.new()
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE),heights=c(4, 4))
## b : mfdfa output
b <- mfdfa
par(mai=rep(0.8, 4))
## 1st plot: Scaling function order Fq (q-order RMS)
p1<-c(1,which(q==0),which(q==q[length(q)]))
plot(log2(scale),log2(b$Fqi[,1]),  pch=16, col=1, axes = F, xlab = "s (days)",
     ylab=expression('log'[2]*'(F'[q]*')'), cex=1, cex.lab=1.6, cex.axis=1.6,
     main= "Fluctuation functionFq",
     ylim=c(min(log2(b$Fqi[,c(p1)])),max(log2(b$Fqi[,c(p1)]))))

lines(log2(scale),b$line[,1], type="l", col=1, lwd=2)
grid(col="midnightblue")
axis(2)
lbl<-scale[c(1,floor(length(scale)/8),floor(length(scale)/4),
             floor(length(scale)/2),length(scale))]
att<-log2(lbl)
axis(1, at=att, labels=lbl)
for (i in 2:3){
  k<-p1[i]
  points(log2(scale), log2(b$Fqi[,k]),  col=i,pch=16)
  lines(log2(scale),b$line[,k], type="l", col=i, lwd=2)
}
legend("bottomright", c(paste('q','=',q[p1] , sep=' ' )),cex=2,lwd=c(2,2,2),
 bty="n", col=1:3)

## 2nd plot: q-order Hurst exponent
plot(q, b$Hq, col=1, axes= F, ylab=expression('h'[q]), pch=16, cex.lab=1.8,
    cex.axis=1.8, main="Hurst exponent", ylim=c(min(b$Hq),max(b$Hq)))
grid(col="midnightblue")
axis(1, cex=4)
axis(2, cex=4)

## 3rd plot: q-order Mass exponent
plot(q, b$tau_q, col=1, axes=F, cex.lab=1.8, cex.axis=1.8,
     main="Mass exponent",
     pch=16,ylab=expression(tau[q]))

grid(col="midnightblue")
axis(1, cex=4)
axis(2, cex=4)

## 4th plot: Multifractal spectrum
plot(b$spec$hq, b$spec$Dq, col=1, axes=F, pch=16, #main="Multifractal spectrum",
     ylab=bquote("f ("~alpha~")"),cex.lab=1.8, cex.axis=1.8,
     xlab=bquote(~alpha))
grid(col="midnightblue")
axis(1, cex=4)
axis(2, cex=4)
x1=b$spec$hq
y1=b$spec$Dq
rr<-poly_fit(x1,y1,4)
mm1<-rr$model1
mm<-rr$polyfit
x2<-seq(0,max(x1)+1,0.01)
curv<-mm[1]*x2^4+mm[2]*x2^3+mm[3]*x2^2+mm[4]*x2+mm[5]
lines(x2,curv, col="red", lwd=2)
reset()
legend("top", legend="MFDFA Plots", bty="n", cex=2)