# How to simulate Multifractional Brownian Motion?

The maximum I would ask, is hints, how to model MFBMotion.

Firstly, i wanted to simulate in MATLAB that code:

nT  = ceil(T/step);
s,c] = cov2corr(sigma);
s   = s * sqrt(step);
sigma = corr2cov(s,c);
Sd  = chol(sigma);
dW  = randn(nT,length(S0)*nb_traj);
Sd  = repmat({Sd},1,nb_traj);
Sdb = blkdiag(Sd{:});
dWc = dW * Sdb;
c   = repmat(1+mu*step, nT, nb_traj);
S   = [repmat(S0, 1, nb_traj); repmat(S0, nT, nb_traj) .* cumprod(c + dWc)];
S   = reshape(S, [nT+1, size(S0,2), nb_traj]);
if nargout > 1
t = [0;step * (1:nT)'];
end


Don't know, how to convert it to Mathematica. Plus what to do with that co-variation function

• And the maximum we'll ask is that you show your Mathematica attempts. Then we can unlock your question. – J. M.'s ennui Oct 23 '15 at 10:17
• Rather than trying to just translate, why not go back to the original sources and develop it from them. It will give you a better understanding of what Mandelbrot does and I'd wager a more elegant Mma implementation. See answer to this question for references. Mandelbrot essentially applies random selection to both time and volatility. Everything else just details. – Jagra Oct 23 '15 at 14:05
• Hmmm...? After reading through link, I don't know whether you want to look at Mandelbrot's implementation o this sort of thing or something else entirely. The paper does reference him. What do you want to accomplish? What problem do you want to solve or model with this? – Jagra Oct 23 '15 at 23:20
• The matlab code is missing a left square bracket on line 2. You need to indicate what are the inputs and supply an example for those values. cov2corr appears to be a local function (doesn't exist in my version of Matlab, 8.1). Similar remarks for corr2cov. Can you add those functions or elaborate if they are new in your version? – Jack LaVigne Oct 25 '15 at 16:17

To help get you started I attempted to translate equation (5) from your snapshot into Mathematica.

Step 1

I(H) from the snapshot becomes:

i[h_] := Piecewise[{
{Gamma[1 - 2 h]/h Sin[π/2 (1 - 2 h)], 0 < h < 1/2},
{Gamma[2 (1 - h)]/(h (2 h - 1)) Sin[π/2 (2 h - 1)],
1/2 < h < 1},
{π, h == 1/2}
}]


This is a function that can be plotted in the known range zero to one.

Plot[i[h], {h, 0, 1}]


Step 2

g(H1,H2) from the snapshot becomes:

g[h1_, h2_] := i[(h1 + h2)/2]/Sqrt[i[h1] i[h2]]


It can be plotted using the known range zero to one. Exclusions -> None was set so there are no interruptions where the derivative is discontinuous.

Plot3D[g[h1, h2], {h1, 0, 1}, {h2, 0, 1}, Exclusions -> None]


Step 3

In creating the covariance between BH1(t1) and BH2(t2) I will assume that t is greater than zero so the absolute times can be removed from all but the difference in time:

covB1B2[h1_, h2_, t1_, t2_] :=
g[h1, h2]/2*(t1^(h1 + h2) + t2^(h1 + h2) - Abs[t1 - t2]^(h1 + h1))


In order to plot this we need to make assumptions about the h1 and h2 values. Further we need to know the range of the time (I have no knowledge of this, nor do I know the time units). I chose zero to one for the time.

The function Manipulate can be used to dynamically set the h1 and h2 parameters.

Manipulate[
Plot3D[covB1B2[h1, h2, t1, t2], {t1, 0, 1}, {t2, 0, 1},
Exclusions -> None],
{{h1, 0.3}, 0, 1, Appearance -> "Open"},
{{h2, 0.7}, 0, 1, Appearance -> "Open"}
]


Change the parameters by using the slider or typing in the values.