A fractional Brownian motion is a Gaussian process satisfying
b) is Gaussian
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.
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
so how to achieve ? The property of circulant matrix shows that
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
r=nan(n+1,1); r(1) = 1;
r(k+1) = 0.5*((k+1)^(2*H) – 2*k^(2*H) + (k-1)^(2*H));
r=[r; r(end-1:-1:2)]; % first rwo of circulant matrix
lambda=real(fft(r))/(2*n); % eigenvalues
W = n^(-H)*cumsum(real(W(1:n+1))); % rescale
%Terminal=W(n+1);%terminal value of fBM