# A simple procedure to simulate multifractional Brownian motion paths

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.

• I agree with you: this sounds very strange to me. Could it be that the authors meant for you to correlate the $k$-th increment of the mBm as if sampled from a fBm with covariance up to time $k$? The method they propose would not really give mBm paths as it does not really respect the mBm covariance structure. May 9, 2022 at 19:33
• @JoseAvilez How can we correlate the $k$-th increment of the mBm in the way you described? May 10, 2022 at 5:24
• @JoseAvilez I think I got it. Every algorithm which generates Brownian motion paths creates a normal random vector for each path. Same for fBm algorithms. Also the function written in the code above, i.e. fbm(hs(k), n), creates internally a normal random vector, call it rand_vec, each time it is called. So I guess that rand_vec must be created before the for loop, and then passed to the fbm function, i.e. fbm(hs(k), n, rand_vec). I tried it and it works. However Peltier and Vehel don't mention this fact in their procedure, maybe because it is obvious? May 11, 2022 at 8:34