Since writing about the simulation of fractional Brownian motion, I’ll spend an hour to write RKHS.

Definition: A Hilbert space $\mathcal{H}$ is a RKHS if $|\mathcal{F}_t[f]|=|f(t)|\leq M\|f\|_{\mathcal{H}}, \forall f\in \mathcal{H}$

Theorem: If $\mathcal{H}$ is a RKHS, then for each $t \in X$ there exists a function $K_t \in \mathcal{H}$ (called the representer of $t$ ) with the reproducing property

$\mathcal{F}_t[f]=<K_t,f>_{\mathcal{H}}=f(t), \forall f\in\mathcal{H}$

Therefore, $K_t(t^{'})=<K_t,K_{t^{'}}>_{\mathcal{H}}$ .

Definition: $K: X\times X\rightarrow\mathbb{R}$ is a reproducing kernel if it’s symmetric and positive definite.

Theorem: A RKHS defines a corresponding reproducing kernel. Conversely, a reproducing kernel defines a unique RKHS.

Once we have the kernels, if $f(\cdot)=\sum \alpha_i K(t_i,\cdot), g(\cdot)=\sum \beta_i K(t^{'}_i,\cdot)$ ,

then $<f,g>_{\mathcal{H}}=\sum\sum \alpha_i \beta_j K(t_i, t^{'}_j)$

When it comes to the fractional Brownian motion

Theorem: For fBm RKHS $K(x,x^{'})=R(x,x^{'})$ is symmetric and positive definite, , there exists $K^H(x,\cdot)$ s.t $K(x,x^{'})=\int K^H(x,y)K^H(x^{'},y)dy$

Take Wiener integral as an example:

$\begin{pmatrix} \textit{dual space}&\textit{inner product}&E&\leftrightarrow&<\cdot,\cdot>_{\mathcal{H}}&\leftrightarrow&<\cdot,\cdot>_{L^{2}}\\ \delta_{t}(\cdot)& &B^H_t&\leftrightarrow&R(\cdot,t)&\leftrightarrow&K^H(t,\cdot)\\ f(\cdot)& &\textit{"}\int_0^T f(t)B^H_tdt\textit{"}&\leftrightarrow&K^H\circ(K^H)^{*}f&\leftrightarrow&(K^H)^{*}f \end{pmatrix}$

where $K^H\circ f(t)=\int_0^T K^H(t,s)f(s)ds$ ,

$(K^H)^{*}\circ f(t)=\int_0^T K^H(s,t)f(s)ds$ and

$K^H\circ(K^H)^{*}f(t)=\int_0^T R(t,s)f(s)ds$

For example, when $H=\frac{1}{2}$ , $K^H(t,s)=1_{[0,t]}(s),$ , we can check

$E(B_tB_s)=<t\wedge\cdot, s\wedge\cdot>_{\mathcal{H}}=<1_{[0,t]}(\cdot),1_{[0,s]}(\cdot)>_{L^2}$

Claim: By previous theorem, the reproducing kernel $R(t,s)$ uniquely defines a RKHS $L(R(t,\cdot))$ with the inner product $\mathcal{H}$

how to simulate a fractional Brownian motion?

Posted on April 7, 2015 by What's your question today? __Mr Curiosity

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