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