# how to simulate a fractional Brownian motion?

A fractional Brownian motion $B^H_t$is 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