Given $(B_H(t))_{t\ge0}$ fractional Brownian motion (fBm) with Hurst parameter $H\in(0,1)$, that is, $B_H(0)=0$, $\mathbb E[B_H(t)]=0$ for all $t$, and for $s,t\geqslant 0$ $B_H$ has the covariance function $$\mathbb E[B_H(t)B_H(s)] = \frac12(|t|^{2H}+|s|^{2H}-|t-s|^{2H}), $$ compute the integral of the fBm w.r.t. time:
$$I = \int_{h}^{2h} B_H(s) ds. $$
Does an exact method to compute $I$ exist?
When dealing with classical Brownian motion $(B(t))_{t\ge0}$, I remember that its integral w.r.t. time could be approximated with the Riemann sum.
Is it possible to do the same with the fBm? I would say yes since $B_H$ is a bounded and continuous Gaussian process.
With the change of variable $x = s-h,\ dx = ds$, we can write $I$ as
$$ I = \int_0^{h} B_H(x+h)dx \overset{?}{\approx} \sum_{i=0}^{N-1} (t_{i+1}-t_i)B_H(t_i+h) $$
where $0 = t_0 < t_1 < ... < t_N = t$ is a partition of $[0,h]$.
Is this approximation correct? Is there any better approximation?