Simulating fractional Brownian motion

Question

I am trying to simulate fBm with its integral representation. I know that there are faster methods out there, but i would like to play around with the kernel function inside the integral. enter link description here My appraoch was to simulate the stochastic integral cummulatively.

set.seed(100)
dt = 0.01

T = seq(-50,10,0.1)

n=length(T)
m=length(T)
Gamma = matrix(0, n, m)

H=0.9
exponent <- function(a, pow) (abs(a)^pow)*sign(a)


for(j in 1:length(T)){
  zeile = numeric(length(T)) #resetting our path for each j
  y = sqrt(dt)*rnorm(n=1, mean = 0, sd=1) #normal distributred r.V.
  zeile[1] = (max(exponent(T[j] - T[1],H-0.5),0) - max(exponent(-T[1],H-0.5),0))*y #first entry of one path
    for(i in 1:(length(T)-1)){
      y1 = sqrt(dt)*rnorm(n=1, mean = 0, sd=1)
      y2 = sqrt(dt)*rnorm(n=1, mean = 0, sd=1)
      zeile[i+1] = zeile[i] + max(exponent(T[j] - T[i],H-0.5),0)*y1 - max(exponent(-T[i],H-0.5),0)*y2
    }
  Gamma[j,] = zeile
}

normalV = rnorm(length(T),mean =0,sd=1)
path = Gamma%*%normalV

plot(T, path, type="l")

T is the interval over which we will plot. The first for loop is there to go over each time point and fix it in the second for loop in which we fill up one row of the matrix. After we have our matrix, I thought I have to multiply it by a N(0,1) vektor in order to get one path of fBm. Clearly I do not.