Gaussian_Process_via_Cholesky.Rmd

---
title: 'Gaussian process simulations with Cholesky method'
author: "Céline Nguyen-Tu & Raphaël Mignot"
date: "12 April 2019"
runtime: shiny
output:
  html_notebook:
    number_sections: no
    toc: yes
  html_document:
    encoding : "utf-8"
    df_print: paged
    toc: yes
  pdf_document:
    toc: yes
linkcolor: red
header-includes:
- \usepackage{bbm}
- \usepackage{amsmath}
- \usepackage{amssymb}
- \usepackage{amsthm}
- \usepackage{natbib}
- \usepackage(mathtools)
---

```{r setup, include=FALSE}
# Required packages
packages <- c("magrittr", "gridExtra", "grid", "latex2exp", "MASS", "knitr", "shiny")
# Install necessary packages if not installed yet
if (any(!requireNamespace(packages, quietly = TRUE))) {
  install.packages(packages, dependencies = TRUE)
}
# Load necessary packages
lapply(packages, require, character.only = TRUE)
```

```{r global_options, echo=FALSE}
knitr::opts_chunk$set(fig.align = 'center',fig.pos = 'h')
library(magrittr)
library(gridExtra)
library(grid)
library(latex2exp)
library(MASS)
library(knitr)
library(shiny)
```

## Gaussian processes 

A vector $\mathbf{\xi} = (X_1 \dots X_n) \in \mathbb{R}^2$ is said to be normal or Gaussian if for all $a_1, \dots , a_n \in \mathbb{R}$, the real random variable $\sum_{i=1}^{n}a_iX_i$ is normally distributed. $\mathbf{\xi}$ can thus be written in the following notation : $\xi \sim \mathcal{N}_n (\textbf{m},\mathbf{\Gamma})$ <br/>
with n-dimentional mean vector $\textbf{m} = \mathbb{E}[\mathbf{\xi}]  = \begin{pmatrix} \mathbb{E}[X_1]\\ \vdots  \\ \mathbb{E}[X_n] \end{pmatrix}$
and $n \times n$ covariance matrix $\mathbf{\Gamma} = \begin{pmatrix} &  & \\  & \text{Cov}(X_i, X_j) & \\  &  & \end{pmatrix}_{i=1 \dots n, j=1 \dots n }$
<br/>
<br/>

A $\mathbb{R}^d$-valued discrete-time stochastic process $\mathbf{X} = \left \{ X_t \; ; \; 0 \leq t < +\infty \right \}$ is called Gaussian if, for any integer $k \geq 1$ and real numbers $0 \leq t_1 < t_2 < \dots < t_k < +\infty$, the random vector $(X_{t_1}, \dots, X_{t_k})$ has a joint Gaussian distribution. In other terms, $\mathbf{X} = (X_t)_t$ is Gaussian if all linear combination $a_1X_{t_1} + \cdots + a_kX_{t_k}$ is normal $\forall k \in \mathbb{N}, \; t_1,\dots, t_k  \in T, \; a_1, \dots, a_k \in \mathbb{R}$. Since we are working on discrete-time processes, $T \in \mathcal{P}(\mathbb{N})$.

The finite-dimentional distributions of a Gaussian process $\mathbf{X}$ are determined by its mean function $m(t) = \mathbb{E}[X_t]$; $t \geq 0$, and its covariance function $K(s,t)= \text{Cov}(X_s, X_t)$; $0 < t, s < +\infty$. We can write $\mathbf{X} \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})$ <br/>
where $\mathbf{\mu} = \mathbb{E}[\mathbf{X}]  = \begin{pmatrix} \mathbb{E}[X_{t_1}]\\ \vdots  \\ \mathbb{E}[X_{t_n}] \end{pmatrix}$ and $\mathbf{\Sigma} = \begin{pmatrix} K(X_{t_1}, X_{t_1})  & \cdots  & K(X_{t_1}, X_{t_n}) \\ \vdots & \ddots & \vdots \\ K(X_{t_n}, X_{t_1})  & \cdots &  K(X_{t_n}, X_{t_n})\end{pmatrix}$

## Presentation of the Cholesky method 

André-Louis Cholesky was a French military officer and mathematician who gave his name to the Cholesky factorization in linear algebra. The latter is a decomposition of a Hermitian positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose. The Cholesky decomposition of the positive-definite matrix $\textbf{A}$ is a decomposition of the form $\textbf{A} = \textbf{LL}^{T}$ where $\textbf{L}$ is a lower triangular matrix with real and positive diagonal entries and $\textbf{L}^{T}$ is the conjugate transpose of $\textbf{L}$. This decomposition is unique. 

Our aim here is to generate a nondegenerate finite Gaussian process $\mathbf{X}$ with mean $\mathbf{\mu}$ and covariance $\mathbf{\Sigma}$ using the Cholesky method through the following steps: <br/>

1. Compute the Cholesky decomposition of $\mathbf{\Sigma}$, i.e. $\mathbf{\Sigma}= \textbf{LL}^{T}$ using the r function `chol`(warning: `chol` returns a upper triangular matrix).
2. Generate a multivariate Gaussian random variable $\mathbf{U} \sim \mathcal{N}_n(0, I_n)$.
3. Then $\mathbf{X} = \mathbf{\mu} \; + \; \mathbf{L}\mathbf{U} \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})$.

We can easily prove that $\mathbf{X}$ has the right distribution: <br/>
- $\mathbf{X} = \mathbf{\mu} \; + \; \mathbf{L}\mathbf{U}$ is Gaussian as it is a linear combination of the Gaussian vector $\mathbf{U}$ <br/>
- $\mathbb{E}[\mathbf{X}]= \mathbb{E}[\mathbf{\mu} \; + \; \mathbf{L}\mathbf{U}] = \mathbb{E}[\mu] + \mathbb{E}[ \mathbf{L}\mathbf{U} ] = \mathbf{\mu} \; + \; \mathbf{L} \cdot \mathbb{E}[\mathbf{U}] =  \mathbf{\mu}$ because $\mathbb{E}[\mathbf{U}] = 0_n$<br/>
- $\mathbb{V}(\mathbf{X})= \mathbb{V}(\mathbf{\mu} \; + \; \mathbf{L}\mathbf{U}) = \mathbf{L} \cdot I_n \cdot \mathbf{L}^{T} = \textbf{LL}^{T} = \mathbf{\Sigma}$ because $\mathbb{V}(\mathbf{U}) = I_n$

The r function `GaussianProcess` simulates a Gaussian process $\mathbf{X} = \left \{ X_t \right \}_{t=t_0}^{t_n}$ from the following arguments: <br/>
* `from`and `to` are numerics for minimum and maximum values of time range <br/>
* `m` is the mean function of the Gaussian process <br/>
* `K` is the covariance function of the Gaussian process <br/>
* `start` is the real value of starting position of the process <br/>
* `n` is the positive integer for the number of points in the process to simulate <br/>
In order to perform the Cholesky method, the program completes a discretization of the chosen time interval. It returns a data frame object of two components: a list with values $t_0, \dots, t_n$ and a list of the corresponding values $X_{t_0}, \dots, X_{t_n}$. If the function is called with its default arguments, it simulates a Brownian motion. 

```{r, include=FALSE}
GaussianProcess = function(from = 0, to = 1, m = function(t) {return(0)}, K = function(s, t) {return(min(s, t))}, start = 0, n = 1000) {
  dt <- (to-from)/n
  tps <- seq(from = from+dt, to = to, length.out = n)
  Sigma <- matrix(nrow=n, ncol=n)
  for (i in 1:n){
    for (j in 1:n){
      Sigma[i,j] <- K(tps[i], tps[j])
    }
  }
  Mu <- sapply(tps, function(t) m(t))
  L <- t(chol(Sigma))
  U <- rnorm(n)
  X <- Mu + L%*%U
  return(data.frame("t" = c(from, tps), "Xt" = c(start, X)))
}
```

## Application : Brownian motion

The most famous Gaussian process is the Brownian motion : it is used to describe a phenomenon in which some quantity is constantly undergoing small, random fluctuations. 
Let $T = \mathbb{R}_+$, a real-valued stochastic process $(W_t)_{t \geq 0}$ is a Brownian motion starting at $x \in \mathbb{R}_+$ if the following conditions are satisfied: <br/>
1. $W_0 = x$ <br/>
2. The function $t \mapsto W_t$ is continuous almost everywhere <br/>
3. For all $0 \leq s \leq t$, $W_t - W_s$ is independent of $W_r$ for all $r \leq s$ <br/>
4. For all $0 \leq s \leq t$, $W_t - W_s$ is normally distributed with with expectation value zero and variance $t-s$ <br/>

The Brownian motion $(W_t)_{t \geq 0}$ is said to be standard if $W_0 = 0$, in which case it is often called the Wiener process. 

All standard Brownian motion $(W_t)_{t \geq 0}$ can be defined by <br/>

<center>
$\mathbb{E}[W_t]=0$  and  $K(s,t) = \text{min}(s,t)$ , $s,t \in T$</center>
<br/>
**Proof** : Let us assume $(W_t)_{t \geq 0}$ is a Brownian motion and let $0 \leq s \leq t$. Then, we have $\mathbb{E}[W_t]= \mathbb{E}[W_t-0] = \mathbb{E}[W_t-W_0] = 0$ because $W_t-W_0 \sim \mathcal{N}(0,t)$. Moreover, using the fact that $W_t - W_s$ is independent of $\sigma (W_r, r\leq s)$, and in particular of $W_s$, then $\mathbb{E}[W_s(W_t-W_s)]=\mathbb{E}[W_s]\mathbb{E}[W_t-W_s]=0$ and so $\text{Cov}(W_s, W_t) = \mathbb{E}[W_sW_t]-\mathbb{E}[W_s]\mathbb{E}[W_t] = \mathbb{E}[W_s(W_t-W_s)+W_s^2]= \mathbb{E}[W_s(W_t-W_s)] + \mathbb{E}[W_s^2]= \mathbb{E}[W_s^2]= s = \text{min}(s,t)$.


```{r, fig.width=5, fig.height=2, include=FALSE}
# Exemples sans library(shiny)
par(mfrow=c(1,2))
BM <- GaussianProcess(from=0, to=10, start=0, n=1000)
plot(BM$t,BM$Xt, type='l', main="Standard Brownian Motion", cex.main=1,
     xlab="t", ylab=bquote(W[t]), col='steelblue')
BM1 <- GaussianProcess(from=0, to=10, start=0, n=1000)
plot(BM1$t,BM1$Xt, type='l', main="Standard Brownian Motion", cex.main=1,
     xlab="t", ylab=bquote(W[t]), col='steelblue', ylim=c(-8,8))
BM2 <- GaussianProcess(from=0, to=10, start=0, n=1000)
points(BM2$t,BM2$Xt, type='l', col='goldenrod')
BM3 <- GaussianProcess(from=0, to=10, start=0, n=1000)
points(BM3$t,BM3$Xt, type='l', col='firebrick')
BM4 <- GaussianProcess(from=0, to=10, start=0, n=1000)
points(BM4$t,BM4$Xt, type='l', col='chocolate1')
BM5 <- GaussianProcess(from=0, to=10, start=0, n=1000)
points(BM5$t,BM5$Xt, type='l', col='chartreuse4')
BM6 <- GaussianProcess(from=0, to=10, start=0, n=1000)
points(BM6$t,BM6$Xt, type='l', col='darkmagenta')
```

```{r motion, echo=FALSE}
ui <- fluidPage(
  hr(),
  fluidRow(
    column(4, numericInput(inputId="Nb", label="Number of processes (up to 6)", value=1, min=1, max=6, step=1)),
    column(4, offset=1.5, sliderInput(inputId="Time", label = "Time range", value=10., min=1., max=100., step=1.)),
    column(4, numericInput(inputId="x", label = "Enter starting position of process", value = 0, min=-100, max=100, step=1))),
    submitButton("Generate"),
    plotOutput("motion"),
  hr()
)

server <- function(input,output){
    output$motion <- renderPlot({
        Time = input$Time
        Nb = input$Nb
        x = input$x
        colors = c('steelblue', 'goldenrod', 'chocolate1', 'chartreuse4', 'darkmagenta')
        K = function(s,t){return(min(s,t))}
        m = function(t){return(x)}
        BB <- GaussianProcess(from=0, to=Time, m=m, K=K, start=x, n=1000)$t
        mini <- 100
        maxi <- -100
        for (k in 1:Nb) {
          new <- GaussianProcess(from=0, to=Time, m=m, K=K, start=x, n=1000)$Xt
          BB <- rbind(BB,new)
          mini <- min(mini, min(new))
          maxi <- max(maxi, max(new))
        }
        par(mar=c(4,4.5,4,4.5))
        plot(BB[1,],BB[2,], type='l', main=paste(c("Brownian motion with starting position x = ", x), collapse=" "), cex.main=1, xlab="t", ylab=bquote(W[t]), col='firebrick', ylim=c(mini, maxi))
        abline(x,0, lty=5, col="grey")
        for (i in 1:Nb-1){
            points(BB[1,],BB[i+2,], type='l', col=colors[i])
          }
        })
}
shinyApp(ui=ui, server=server, options = list(width = "100%", height = 625))
```

## Application : Brownian bridge

Let $T=[0,\alpha] \subset \mathbb{R}_+$, a Brownian bridge is a Wiener process $(W_t)_{t \in T}$ defined on the interval $T$ and conditioned on the event $W_\alpha=0$. We can write this process $(W_t \mid W_\alpha=0)_{t\in [0,\alpha]}$ and call it the tied-down Wiener process. If $\alpha = 1$, then it is commonly called the standard Brownian bridge. The Brownian bridge is an example of a Gaussian process whose increments are not independent.
The standard Brownian bridge $(Z_t)_{t\in [0,1]}$ is defined by $Z_t = W_t-tW_1$ for all $t \in [0,1]$ where $(W_t)_{t \in [0,1]}$ is a Wiener process. <br/>

All standard Brownian bridge can also be defined by <br/>

<center>
$\mathbb{E}[Z_t]=0$ and $K(s,t) = \text{min}(s,t) -st$ , $s,t \in T$ </center>
<br/>
**Proof** : Let us assume $(Z_t)_{t\in [0,1]}$ is a standard Brownian bridge and let $0 \leq s \leq t \leq 1$. Thus, we have $\mathbb{E}[Z_t]= \mathbb{E}[W_t-tW_1] =\mathbb{E}[W_t]-t\mathbb{E}[W_1] = 0$ because $\mathbb{E}[W_t]=0$ for all $t\geq 0$. Moreover, $\text{Cov}(Z_s, Z_t) = \mathbb{E}[(Z_s-\mathbb{E}[Z_s])(Z_t-\mathbb{E}[Z_t])]= \mathbb{E}[(W_s-sW_1)(W_t-tW_1)] = \mathbb{E}[W_sW_t-tW_sW_1 - sW_tW_1 + stW_1^2]$
$\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; = \text{min}(s,t) - t \; \text{min}(s,1) - s \; \text{min}(t,1) + st = \text{min}(s,t) -st-st+st = \text{min}(s,t)-st$

```{r, include=FALSE}
K = function(s,t){return(min(s,t)-s*t)}
```

```{r, fig.width=5, fig.height=2, include=FALSE}
# Exemples sans library(shiny)
par(mfrow=c(1,2))
BM <- GaussianProcess(from=0, to=0.999, K=K, start=0, n=1000)
plot(BM$t,BM$Xt, type='l', main="Brownian Bridge", cex.main=1,
     xlab="t", ylab=bquote(W[t]), col='steelblue')
BM1 <- GaussianProcess(from=0, to=0.999, K=K, start=0, n=1000)
plot(BM1$t,BM1$Xt, type='l', main="Brownian Bridge", cex.main=1,
     xlab="t", ylab=bquote(W[t]), col='steelblue', ylim=c(-1.5,1.5))
BM2 <- GaussianProcess(from=0, to=0.999, K=K,start=0, n=1000)
points(BM2$t,BM2$Xt, type='l', col='goldenrod')
BM3 <- GaussianProcess(from=0, to=0.999, K=K, start=0, n=1000)
points(BM3$t,BM3$Xt, type='l', col='firebrick')
BM4 <- GaussianProcess(from=0, to=0.999, K=K, start=0, n=1000)
points(BM4$t,BM4$Xt, type='l', col='chocolate1')
BM5 <- GaussianProcess(from=0, to=0.999, K=K, start=0, n=1000)
points(BM5$t,BM5$Xt, type='l', col='chartreuse4')
BM6 <- GaussianProcess(from=0, to=0.999, K=K, start=0, n=1000)
points(BM6$t,BM6$Xt, type='l', col='darkmagenta')
```

```{r bridge, echo=FALSE}
ui <- fluidPage(
  hr(),
  fluidRow(
    column(4, numericInput(inputId="Nb", label="Number of processes (up to 6)", value=1, min=1, max=6, step=1)),
    column(4, offset=1.5, sliderInput(inputId="Time", label = "Time range", value=1., min=1., max=100., step=1.)),
    column(2, submitButton("Generate"))),
    plotOutput("bridge"),
  hr()
)

server <- function(input,output){
    output$bridge <- renderPlot({
        Time = input$Time
        Nb = input$Nb
        colors = c('steelblue', 'goldenrod', 'chocolate1', 'chartreuse4', 'darkmagenta')
        K = function(s,t){return(min(s,t)-(s*t)/Time)}
        BB <- GaussianProcess(from=0, to=Time-0.001, K=K, start=0, n=1000)$t
        mini <- 0
        maxi <- 0
        for (k in 1:Nb) {
          new <- GaussianProcess(from=0, to=Time-0.001, K=K, start=0, n=1000)$Xt
          BB <- rbind(BB,new)
          mini <- min(mini, min(new))
          maxi <- max(maxi, max(new))
        }
        par(mar=c(4,4.5,4,4.5))
        plot(BB[1,],BB[2,], type='l', main=paste(c("Brownian Bridge for T = [ 0 ,", Time, "]"), collapse=" "), cex.main=1, xlab="t", ylab=bquote(Z[t]), col='firebrick', ylim=c(mini, maxi))
        abline(0,0, lty=5, col="grey")
        for (i in 1:Nb-1){
            points(BB[1,],BB[i+2,], type='l', col=colors[i])
          }
        })
}
shinyApp(ui=ui, server=server, options = list(width = "100%", height = 590))
```


## Application : Fractional Brownian motion

Let $T = \mathbb{R}_+$, the fractional Brownian motion $(B^H_t)_{t \geq 0}$, also called fractal Brownian motion, is a generalization of Brownian motion. It is an example of a Gaussian process whose increments do not need to be independent. The fractional Brownian motion is defined as the only Gaussian process which starts at zero and such that <br/>

<center>
$\mathbb{E}[B^H_t]=0$ and $K(s,t) = \frac{1}{2}(|s|^{2H}+|t|^{2H} - |s-t|^{2H})$ </center>

where $H \in ]0,1[$ is a real constant value called the Hurst index. <br/>

* If $H=\frac{1}{2}$, the fractional Brownian motion is a standard Brownian motion. <br/>
* If $H > \frac{1}{2}$, this process corresponds to a fractional integral of the Brownian motion. <br/>
* If $H < \frac{1}{2}$, it corresponds to a fractional derivative of the Brownian motion.

```{r, include=FALSE}
H=0.8
K = function(s,t){return(0.5*(abs(s)^(2*H)+abs(t)^(2*H)-abs(s-t)^(2*H)))}
```

```{r, fig.width=5, fig.height=5, include=FALSE}
# Exemples sans library(shiny)

par(mfrow=c(2,2))

fBM <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
plot(fBM$t,fBM$Xt, type='l', main="Fractional Brownian Motion with H=0.8", cex.main=1, xlab="t", ylab=bquote(W[t]), col='steelblue')

fBM1 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
plot(fBM1$t,fBM1$Xt, type='l', main="Fractional Brownian Motion with H=0.8", cex.main=1,xlab="t", ylab=bquote(W[t]), col='steelblue', ylim=c(-10,10))
fBM2 <- GaussianProcess(from=0, to=10, K=K,start=0, n=1000)
points(fBM2$t,fBM2$Xt, type='l', col='goldenrod')
fBM3 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM3$t,fBM3$Xt, type='l', col='firebrick')
fBM4 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM4$t,fBM4$Xt, type='l', col='chocolate1')
fBM5 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM5$t,fBM5$Xt, type='l', col='chartreuse4')
fBM6 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM6$t,fBM6$Xt, type='l', col='darkmagenta')

H=0.1
K = function(s,t){return(0.5*(abs(s)^(2*H)+abs(t)^(2*H)-abs(s-t)^(2*H)))}
fBM1 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
plot(fBM1$t,fBM1$Xt, type='l', main="Fractional Brownian Motion with H=0.1", cex.main=1,xlab="t", ylab=bquote(W[t]), col='steelblue', ylim=c(-10,10))
fBM2 <- GaussianProcess(from=0, to=10, K=K,start=0, n=1000)
points(fBM2$t,fBM2$Xt, type='l', col='goldenrod')
fBM3 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM3$t,fBM3$Xt, type='l', col='firebrick')
fBM4 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM4$t,fBM4$Xt, type='l', col='chocolate1')
fBM5 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM5$t,fBM5$Xt, type='l', col='chartreuse4')
fBM6 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM6$t,fBM6$Xt, type='l', col='darkmagenta')

H=0.4
K = function(s,t){return(0.5*(abs(s)^(2*H)+abs(t)^(2*H)-abs(s-t)^(2*H)))}
fBM1 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
plot(fBM1$t,fBM1$Xt, type='l', main="Fractional Brownian Motion with H=0.4", cex.main=1,xlab="t", ylab=bquote(W[t]), col='steelblue', ylim=c(-10,10))
fBM2 <- GaussianProcess(from=0, to=10, K=K,start=0, n=1000)
points(fBM2$t,fBM2$Xt, type='l', col='goldenrod')
fBM3 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM3$t,fBM3$Xt, type='l', col='firebrick')
fBM4 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM4$t,fBM4$Xt, type='l', col='chocolate1')
fBM5 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM5$t,fBM5$Xt, type='l', col='chartreuse4')
fBM6 <- GaussianProcess(from=0, to=10, K=K, start=0, n=1000)
points(fBM6$t,fBM6$Xt, type='l', col='darkmagenta')
```

```{r frac, echo=FALSE}
ui <- fluidPage(
  hr(),
  fluidRow(
    column(4, numericInput(inputId="Nb", label="Number of processes (up to 6)", value=1, min=1, max=6, step=1)),
    column(4, offset=1.5, sliderInput(inputId="Time", label = "Time range", value=10., min=1., max=100., step=1.)),
    column(4, sliderInput(inputId="H", label = "Choose the value of Hurst constant H", value = 0.8, min=0.01, max=0.99, step=0.01))),
    submitButton("Generate"),
    plotOutput("frac"),
  hr()
)

server <- function(input,output){
    output$frac <- renderPlot({
        H = input$H
        Time = input$Time
        Nb = input$Nb
        colors = c('steelblue', 'goldenrod', 'chocolate1', 'chartreuse4', 'darkmagenta')
        K = function(s,t){return(0.5*(abs(s)^(2*H)+abs(t)^(2*H)-abs(s-t)^(2*H)))}
        fBM <- GaussianProcess(from=0, to=Time, K=K, start=0, n=1000)$t
        mini <- 0
        maxi <- 0
        for (k in 1:Nb) {
          new <- GaussianProcess(from=0, to=Time, K=K, start=0, n=1000)$Xt
          fBM <- rbind(fBM,new)
          mini <- min(mini, min(new))
          maxi <- max(maxi, max(new))
        }
        par(mar=c(4,4.5,4,4.5))
        plot(fBM[1,],fBM[2,], type='l', main=paste(c("Fractional Brownian Motion with H =", H), collapse=" "), cex.main=1, xlab="t", ylab=bquote(B[t]^H), col='firebrick', ylim=c(mini, maxi))
        abline(0,0, lty=5, col="grey")
        for (i in 1:Nb-1){
            points(fBM[1,],fBM[i+2,], type='l', col=colors[i])
          }
        })
}
shinyApp(ui=ui, server=server, options = list(width = "100%", height = 625))
```

##### References

[1] Vincent Lemaire & Tabea Rebafka (2019), *Probabilités numériques et Statistiques Computationnelles*, Master 1 Mathematics and Applications, Sorbonne University <br/>
[2] Ioannis Karatzas & Steven E. Shreve (1991), *Brownian Motion and Stochastic Calculus*, Graduate Texts in Mathematics, Springer <br/>
[3] Léonard Gallardo (2008), *Mouvement Brownien et Calcul d'Itô*, Collections Méthodes, Hermann <br/>
[4] Philippe Bougerol (2015), *Modèles stochastique et Applications à la Finance*, Master 1 Mathematics and Applications, Sorbonne University <br/>