# Reproducing Kernel Hilbert Space

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]=_{\mathcal{H}}=f(t), \forall f\in\mathcal{H}$

Therefore, $K_t(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 $_{\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)=_{\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?

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

# comparison

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 ?