A simple procedure to simulate multifractional Brownian motion paths

Question

In a paper by Peltier and Vehel the multifractional Brownian motion (mBm) was defined for the first time, and they also give a procedure to simulate mBm sample paths. Briefly, mBm generalizes the fractional Brownian motion (fBm) which in turn generalizes the Brownian motion (aka Wiener process). fBm generalizes the BM in the sense that it allows the roughness of the paths to be tuned with a constant parameter $H\in(0,1)$ (the higher $H$ is the smoother the curve will be, in particular if $H=0.5$ then fBm = BM). mBm generalizes the fBm in the sense that it allows $H$ to vary along the trajectory, that is $H$ becomes a function of time.

The simulation procedure is described at page 30 of the pdf. For the sake of brevity let $k=i/N$ where $1\le i\le N$ and $N$ is the size of the desired sample, then:

We assume given a function $t \mapsto H(t)$. For each value of $H(k)$, an fBm $B_{H(k)}$ of exponent $H(k)$ is generated. The mBm $W_H$ is then obtained by setting $$W_H(k) = B_{H_k}(k).$$ In practice, any method may be used for generating each fBm.

Assuming that $B_{H(k)}$ and $B_{H_k}$ represent the same thing, what I understand from the procedure above is the following: to generate a mBm sample path of length $N$, we have to generate $N$ fBm sample paths of length $N$, one for each "time" $k$, and each of this fBm is generated using a different Hurst exponent given by $H(k)$. But for each path $B_{H(k)}$ we only need its $k$-th value, i.e. $B_{H(k)}(k)$, while the others $N-1$ values are useless. And this $k$-th value of the $k$-th fBm path will be the $k$-th value of the mBm path. However I think that this reasoning is not correct, because I tried to apply it and the resulting trajectory is not a mBm path (see explanation below the code).

This is the Matlab code

% lenght of the sample path
n = 2^9;
% Hurst function
H = @(t) 0.3 + 0.5*t;
% time steps, i.e. discretization of time interval [0,1] in n equidistant points
ts = linspace(0, 1, n);
% Hurst steps, i.e. discretization of H
hs = H(ts);
% initialize the array
mbm_path = zeros(n,1);
for k = 2:n
    % fBm path generated with the Hurst exponent hs(k) using the Wood-Chan method
    fbm_path = fbm(hs(k), n);
    % the k-th value of the mBm path is given by k-th value of the k-th fBm path
    mbm_path(k) = fbm_path(k);
end

This is a mBm path generated with $H(t) = 0.3 + 0.5t$ for $t\in[0,1]$, but using a different procedure than the one described by Peltier and Levy. (blue line is the mBm path, red line is $H$)

This is a path generated with the same $H$ as before, but using the code above, and it is easy to see that it isn't a mBm path, in fact the roughness of the path increases as $H$ increases (it should be the opposite).

The problem is due to the fact that statistically the distance between two fBm paths will increase as time goes on. Thus at the line mbm_path(k) = fbm_path(k), as k approaches n, it is likely that fbm_path(k) and fbm_path(k+1) are far apart. This effect is evident in particular for $H>0.5$, since in this case the increments of the process are positively correlated hence it is more likely that the trend is preserved. For example the following figure shows two fBm paths with $H=0.8$.

The description of the procedure by Peltier and Vehel is very simple, yet I don't get where I am doing wrong.

Answer 1

The fractional Brownian motions you generate must not be independent. They are generated from the same Gaussian field $(B_H(t))_{(H,t) \in [a,b] \times \mathbb{R}}$ and there exists a dependence according to the values of $H$.