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


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s