require('MASS')

# Discretization 
N <- 500

# Hurstparameter
H <- 0.1

T <- seq(0,N,1)
I <- T/N

mu <- rep(0, N+1)
C <- matrix(nrow=N+1, ncol=N+1)

# Covariance matrix
for(i in 0:N){
  t <- i/N
  for(j in 0:N){
    s <- j/N
    C[i+1,j+1] <- 0.5*(t^(2*H) + s^(2*H) - abs(t-s)^(2*H))
  }
}

X = mvrnorm(n=1, mu, C)
plot(I,X,type='l')