how to simulate a fractional Brownian motion?

A fractional Brownian motion B^H_tis a Gaussian process satisfying

a) B^H_0=0

b) B^H_t is Gaussian

c) E(B^H_tB^H_s)=R(t,s)=\frac{1}{2}(t^{2H}+s^{2H}-|t-s|^{2H})

As for simulation, we just need to simulate the correlated Gaussian (\Delta B^H_{t_1},\dots,\Delta B^H_{t_n})\sim N(0,\Sigma), where

t_i=i*T/n=i \Delta t,

\Delta B^H_{t_i}=B^H_{t_i}-B^H_{t_{i-1}}, i=1,\dots,n,

\Sigma_{i,j}=E(B^H_{t_i}B^H_{t_j})=\frac{\Delta t^{2H}}{2}((i-j+1)^{2H}+(i-j-1)^{2H}-2(i-j)^{2H})

Now the question is just simulating L'\cdot N(0,I_{n\times n}) where L is the Cholesky decomposition of \Sigma, i.e (\Sigma=L'L)

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

Method 2:

We denote \lambda(k)=\frac{(\Delta t)^{2H}}{2}((k+1)^{2H}+(k-1)^{2H}-2k^{2H})=\Sigma_{i,i+k}

Denote a circulant matrix

C=\begin{pmatrix}\lambda(0)&\lambda(1)&\lambda(2)&\dots&\lambda(n-1)&\lambda(n-2)&\dots&\lambda(1)\\ \lambda(1)&\lambda(0)&\lambda(1)&\dots&\lambda(n-2)&\lambda(n-1)&\dots&\lambda(2)\\ \lambda(2)&\lambda(1)&\lambda(0)&\dots&\lambda(n-3)&\lambda(n-2)&\dots&\lambda(3)\\ \vdots&\vdots&\vdots&\ddots&\dots&\dots&\dots&\dots\\    \lambda(n-1)&\lambda(n-2)&\lambda(n-3)&\dots&\lambda(0)&\lambda(1)&\dots&\lambda(n-2)\\    \dots    \end{pmatrix}_{2n-2\times 2n-2}

of which \Sigma is the left up n\times n block matrix.

The idea is, suppose C=DD^*=(D_1+iD_2)*(D_1^T-i D_2^T)where D_1, D_2 \in \mathbb{R}^{2n-2\times 2n-2}, now if we generate i.i.d \quad \epsilon^1,\epsilon^2\sim N(0,I_{2n-2,2n-2}),

We can conclude that

denote X_1+iX_2=D*(\epsilon^1+i\epsilon^2)=D_1\epsilon^1-D_2\epsilon^2+i(D_2\epsilon^1+D_1\epsilon^2)\in \mathbb{R}^{2n-2,1}, then X_1[1:n]\sim N(0,\Sigma) and X_2[1:n]\sim N(0,\Sigma)

let Q is the unitary matrix defined by

Q_{j,k}=(2n-2)^{-1/2}exp(-2i\pi\frac{jk}{2n-2}) and thus D=Q*\Lambda^{\frac{1}{2}}

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

\lambda_k=\sum_{j=1}^{2n-2}c_j exp(-2i\frac{jk}{2n-2}), i.e

\Lambda=Q*\begin{pmatrix}\lambda(0)&\lambda(1)&\lambda(2)&\dots&\lambda(n-1)&\lambda(n-2)&\dots&\lambda(1)\\    \end{pmatrix}^{'} that can be achieved by FFT.

Next, according to the idea stated before, FFT(diag(\Lambda^{\frac{1}{2}})*(\epsilon^1+i\epsilon^2))=Q*diag(\Lambda^{\frac{1}{2}})*(\epsilon^1+i\epsilon^2)=X_1+iX_2

And the first n component of X_1 will be N(0,\Sigma)

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));
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




I was trying to simulate something related to two dimensional brownian motion or fractional brownian motion.

As shown in the graph, can you tell the difference of these two sets of pictures?
At the first glance, they look the same, but actually, their underlying property are very
different. What is it ?