  We introduce a simulation scheme for Brownian semistationary processes, which
is based on discretizing the stochastic integral representation of the process
in the time domain. We assume that the kernel function of the process is
regularly varying at zero. The novel feature of the scheme is to approximate
the kernel function by a power function near zero and by a step function
elsewhere. The resulting approximation of the process is a combination of
Wiener integrals of the power function and a Riemann sum, which is why we call
this method a hybrid scheme. Our main theoretical result describes the
asymptotics of the mean square error of the hybrid scheme and we observe that
the scheme leads to a substantial improvement of accuracy compared to the
ordinary forward Riemann-sum scheme, while having the same computational
complexity. We exemplify the use of the hybrid scheme by two numerical
experiments, where we examine the finite-sample properties of an estimator of
the roughness parameter of a Brownian semistationary process and study Monte
Carlo option pricing in the rough Bergomi model of Bayer et al. [Quant. Finance
16(6), 887-904, 2016], respectively.
