Simulation of brownian motion and fractional brownian motion It's easy to simulate a path of a brownian motion with the method explained in 
Wiener process as a limit of random walk:
import numpy as np
import matplotlib.pyplot as plt
X = 2 * np.random.binomial(1, 0.5, 2*1000*1000) - 1
cumsumX = np.cumsum(X)
n = 1000*1000
x = np.linspace(0, 1, num=n)
Y = 1/np.sqrt(n) * np.array([cumsumX[int(n*t)] for t in x])
plt.plot(x,Y)
plt.show()


Is it possible to simulate a fractional brownian motion with such a  simple method, based on a random walk? (none of the methods explained in this article seem to be random-walk based)
 A: You can, but it is not so easy since the increments of a fractional brownian motion are dependent.
In fact, it is well known that
\begin{eqnarray}\label{Bht}
B^{H}_{t}=\int^{t}_{0}z_{H}(t,s)dB_{s}\,,
\end{eqnarray}where $B_{t}$ is a standard Brownian motion and the deterministic
kernel is given by
\begin{eqnarray}\label{zh}
z_{H}(t,s)=c_{H}\left( \left(\frac{t}{s}\right)^{H-\frac{1}{2}}\left(t-s\right)^{H-\frac{1}{2}}-  \left(H-\frac{1}{2}\right)s^{\frac{1}{2}-H}\int_{s}^{t}u^{H-\frac{3}{2}}
\left(u-s\right)^{H-\frac{1}{2}}du\right){\bf 1}_{(0,t)}(s),
\end{eqnarray}
with the constant
\begin{eqnarray}\label{ch}
c_{H}^2 =\frac{2H\Gamma(\frac{3}{2}-H)}{\Gamma(H+\frac{1}{2})\Gamma(2-2H)}\;,
\end{eqnarray}
and ${\bf 1}_{(0,t)}(s)$  is an
indicator function that takes the value 1 iff $s\in(0,t)$ and is zero otherwise.
Then inspired by Sottinen (2001, FS jounral), for $H\in(0,1)$, we define:
\begin{equation}\label{Bthn}
B_{t}^{H,n}=\sum\nolimits_{i=1}^{\lfloor nt\rfloor} \sqrt{n}
\left(\int_{\frac{i-1}{n}}^{\frac{i}{n}} z_{H}\left(\frac{\lfloor
nt\rfloor}{n},s\right)ds\right)\xi_{i} \;,
\end{equation}
where $\xi_{1},\ldots,\xi_{n}$ is an $n$-tuple of independent symmetric Bernoulli random variables with $\mathbb{P}(\xi_{i}=\pm 1)=\frac{1}{2}$, living on a probability space $(\Omega_{n}, \mathcal {F}_{n},
\mathbb{P}_{n})$.
Now you can similate a fractional brownian motion using $B_{t}^{H,n}$.
