A fractional Brownian motion is a Gaussian process satisfying

a)

b) is Gaussian

c)

As for simulation, we just need to simulate the correlated Gaussian , where

,

,

Now the question is just simulating where is the Cholesky decomposition of , i.e ()

Method 1: find first and then use Cholesky decomposition to compute . This works theoretically, but not practically because it’s too time and space costy to record , not to mention decomposition step.

Method 2:

We denote

Denote a circulant matrix

of which is the left up block matrix.

**The idea is**, suppose where , now if we generate ,

We can conclude that

denote , then and

let Q is the unitary matrix defined by

and thus

so how to achieve ? The property of circulant matrix shows that

, i.e

that can be achieved by FFT.

Next, according to the idea stated before,

And the first n component of will be

Reference : Coeurjolly, JF. Simulation and identification of the fractional Brownian motion

Matlab code:

function []=fBM(n,H)

r=nan(n+1,1); r(1) = 1;

for k=1:n

r(k+1) = 0.5*((k+1)^(2*H) – 2*k^(2*H) + (k-1)^(2*H));

end

r=[r; r(end-1:-1:2)]; % first rwo of circulant matrix

lambda=real(fft(r))/(2*n); % eigenvalues

W=fft(sqrt(lambda).*complex(randn(2*n,1),randn(2*n,1)));

W = n^(-H)*cumsum(real(W(1:n+1))); % rescale

plot((0:n)/n,W);

%Terminal=W(n+1);%terminal value of fBM

end