rm(list=ls())

library(ggplot2)

#generate Brownian motion
t <- 1:1e3
sig2 <- 0.01
x <- rnorm(n = length(t) - 1, sd = sqrt(sig2))
data <- c(0, cumsum(x))

#function to find Hurst exponent
simpleHurst <- function(y){
  sd.y <- sd(y)
  m <- mean(y)
  y <- y - m
  max.y <- max(cumsum(y))
  min.y <- min(cumsum(y))
  RS <- (max.y - min.y)/sd.y
  return(RS)
}

#find Hurst exponent of Brownian motion
df.rs <- data.frame()
for(i in c(1, 2, 4, 8, 16, 32)){
  RS <- 0
  n <- floor((length(data)/i))
  for(j in 0:(i-1)){
    X <- data[(1:n) + j*n]
    RS <- RS + simpleHurst(y = X)
  }

  df.rs <- rbind(df.rs,
                 data.frame("log.t" = log(length(X)),
                            "log.RS"= log(RS / n)))
}

ggplot(data=df.rs,  aes(x=log.t, y=log.RS)) + geom_point(size=I(2))