OUcode.Rmd

---
title: "OUcode"
output: html_notebook
---
# Data Generation
```{r}
## In reality data is given, here we generate it ourself for simulation purpose.
## Note that with the change of measure, the distribution of data Y follows that of a standard Brownian motion.
T = 100
## Lmax decide the finest level of our observations, 2^(-Lmax) is the time interval between each observations
Lmax = 10
Ydis = vector(mode = "numeric",length = T*2^Lmax+1)
bm = vector(mode = "numeric",length = T*2^Lmax)
for(t in 1:(T*2^Lmax+1)){
  bm[t] = rnorm(1, mean = 0, sd = 2^(-Lmax))
  Ydis[t+1] = sum(bm[1:t])
}
```





## True value computation
```{r}
## Use PF at l=10 with 200 replications, to obtain the true value.
## Need to import package e1071 in R

ltemp = 10
# Number of particles N
N = 100 * 2^ltemp
# Time horizon, number of time units
T = 100
# Replications number of PF algorithm to get it's approximate variance, bias and MSE
R = 200
# l here is the discretization level
l = ltemp
# Initialization
Xun = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
ko = vector(mode = "numeric",length = N)
Xun1 = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
estimate_PF = vector(mode = "numeric",length = R)
## replication of PF algorithm
for(h in 1:R){
ess = vector(mode = "numeric",length = T)
Xun[1,]=ko
Xun1[1,]=ko
W = vector(mode = "numeric",length = N)
G = vector(mode = "numeric",length = N)
W = W+1
for (p in 1:T) {
  ## Below is Euler update taking place in time p-1 to p, Xun[p,] represent singals at time p-1 , firstly we use Euler discretization to update Xun[p+1,] from Xun[p,], the goal is to get Xun1[p+1,], which we wiil obtain by perform resample step on Xun[p+1,]
  for(i in 1:N){
  ## particle update using Euler
  xdis = vector(mode = "numeric", length =(2^l+1))
  xdis[1] = Xun1[p,i]
  for (j in 1:(2^l)) {
      dw = rnorm(1, mean=0, sd = sqrt(0.5^l))
      xdis[j+1] = (-xdis[j])*0.5^l + 0.5*dw + xdis[j]
  }
  Xun[p+1,i]=xdis[(2^l+1)]
  ## G weights Function update
  xdd = vector(mode="numeric",length =(2^l))
  for(id in 1:(2^l)){
    xdd[id] = xdis[id]*(Ydis[(p-1)*2^Lmax+id*2^(Lmax-l)+1]-Ydis[(p-1)*2^Lmax+(id-1)*2^(Lmax-l)+1]) - 0.5^(l+1)*(xdis[id]^2)
  }
  G[i]= sum(xdd)
  }
  ## Normalized weight
  W = exp(G-max(G))
  w = W/sum(W)
  ## Calculating Effective Sample Size
  ess[p]=(sum(w^2))^(-1)
  ## If ESS is smaller or equal than the threshold, we perform resample
  if(ess[p]<=(N/2)){
  xindex = vector(mode = "numeric",length = N)
  xindex = rdiscrete(N,W)
  for(i in 1:N){
    Xun1[p+1,i]=Xun[p+1,xindex[i]]
  }
  }
  ## If ESS is larger than the threshold, we don't resample
  if(ess[p]>(N/2)){
    Xun1[p+1,]=Xun[p+1,]
  }
  
}
#ESS = rbind(ESS,ess/N)
estimate_PF[h] = sum(Xun[T+1,]*w)
}
tv = mean(estimate_PF)



```



# OU PF
```{r}
## Gradient Estimator, dYt = Xt*dt + dBt, dXt = (-Xt)*dt + 0.5*dWt (Using change of measure)
date()
estPF = vector(mode = "numeric",length = 8)
MSE_PF = vector(mode = "numeric",length = 8)
var_PF = vector(mode = "numeric",length = 8)
cost_PF = vector(mode = "numeric",length = 8)
for(ltemp in 1:8){
# Number of particles N
N = 100 * 2^ltemp
# Time horizon, number of time units
T = 100
# Replications number of PF algorithm to get it's approximate variance, bias and MSE
R = 200
# l here is the discretization level
l = ltemp
# Initialization
Xun = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
ko = vector(mode = "numeric",length = N)
Xun1 = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
estimate_PF = vector(mode = "numeric",length = R)
## replication of PF algorithm
for(h in 1:R){
ess = vector(mode = "numeric",length = T)
Xun[1,]=ko
Xun1[1,]=ko
W = vector(mode = "numeric",length = N)
G = vector(mode = "numeric",length = N)
W = W+1
for (p in 1:T) {
  ## Below is Euler update taking place in time p-1 to p, Xun[p,] represent singals at time p-1 , firstly we use Euler discretization to update Xun[p+1,] from Xun[p,], the goal is to get Xun1[p+1,], which we wiil obtain by perform resample step on Xun[p+1,]
  for(i in 1:N){
  ## particle update using Euler
  xdis = vector(mode = "numeric", length =(2^l+1))
  xdis[1] = Xun1[p,i]
  for (j in 1:(2^l)) {
      dw = rnorm(1, mean=0, sd = sqrt(0.5^l))
      xdis[j+1] = (-xdis[j])*0.5^l + 0.5*dw + xdis[j]
  }
  Xun[p+1,i]=xdis[(2^l+1)]
  ## G weights Function update
  xdd = vector(mode="numeric",length =(2^l))
  for(id in 1:(2^l)){
    xdd[id] = xdis[id]*(Ydis[(p-1)*2^Lmax+id*2^(Lmax-l)+1]-Ydis[(p-1)*2^Lmax+(id-1)*2^(Lmax-l)+1]) - 0.5^(l+1)*(xdis[id]^2)
  }
  G[i]= sum(xdd)
  }
  ## Normalized weight
  W = exp(G-max(G))
  w = W/sum(W)
  ## Calculating Effective Sample Size
  ess[p]=(sum(w^2))^(-1)
  ## If ESS is smaller or equal than the threshold, we perform resample
  if(ess[p]<=(N/2)){
  xindex = vector(mode = "numeric",length = N)
  xindex = rdiscrete(N,W)
  for(i in 1:N){
    Xun1[p+1,i]=Xun[p+1,xindex[i]]
  }
  }
  ## If ESS is larger than the threshold, we don't resample
  if(ess[p]>(N/2)){
    Xun1[p+1,]=Xun[p+1,]
  }
  
}
#ESS = rbind(ESS,ess/N)
estimate_PF[h] = sum(Xun[T+1,]*w)
}
estPF[l] = mean(estimate_PF)
MSE_PF[l] = (mean(estimate_PF)-tv)^2 + var(estimate_PF)
var_PF[l] = var(estimate_PF)
cost_PF[l] = N*2^l*T
}
date()

# rate for PF here is -2.724165
```


# OU MLPF
```{r}
# Gradient Estimator, dYt = Xt*dt+dBt, dXt = (-Xt)*dt + 0.5*dWt 
# Multilevel Particle Filter structure 
date()
mean_PF = vector(mode = "numeric",length = 8)
MSE_MLPF = vector(mode = "numeric",length = 8)
var_MLPF = vector(mode = "numeric",length = 8)
work_MLPF = vector(mode = "numeric",length = 8)

for(ltemp in 1:8){
# Time horizon, number of time units
T = 100
# l here is the discretization level
L = ltemp
# Replications number of PF algorithm to get it's approximate variance, bias and MSE
R=200
# This vector is used to store the estimate values over 200 replications for MLPF maximal resampling algorithm
estMLPF = vector(mode = "numeric",length = R)

for(h in 1:R){
Aa = vector(mode = "numeric", length = L+1)
for(l in 1:L){
# Initialization
N = 10*ceiling(2^(L-l)*L)
ko = vector(mode = "numeric",length = N)
Xun = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xun1 = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xunc = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xunc1 = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xun[1,]=ko
Xun1[1,]=ko
Xunc[1,]=ko
Xunc1[1,]=ko
wf = vector(mode = "numeric", length = N)
wc = vector(mode = "numeric", length = N)
Wf = vector(mode = "numeric", length = N)
Wc = vector(mode = "numeric", length = N)
Gf = vector(mode = "numeric",length = N)
Gc = vector(mode = "numeric",length = N)
Bb = vector(mode = "numeric", length = N)
# Loop over time
for (p in 1:T) {
  # loop over particles
  for(i in 1:N){
  ## Particle update for one time unit
  dw = vector(mode = "numeric", length =(2^l))
  ## Finer path Euler update
  fdis = vector(mode = "numeric", length =(2^l+1))
  fdis[1] = Xun1[p,i]
  for (j in 1:(2^l)) {
      dw[j] = rnorm(1, mean=0, sd = sqrt(0.5^l))
      fdis[j+1] = (-fdis[j])*0.5^l + 0.5*dw[j] + fdis[j]
  }
  Xun[p+1,i]=fdis[(2^l+1)]
  ##Coarser path Euler update
  cdis = vector(mode = "numeric", length =(2^(l-1)+1))
  cdis[1] = Xunc1[p,i]
  for(j in 1:(2^(l-1))){
    cdis[j+1] = (-cdis[j])*0.5^(l-1) + 0.5*(dw[2*j]+dw[2*j-1]) + cdis[j]
  }
  Xunc[p+1,i]=cdis[(2^(l-1)+1)]
  ## G Function calculation
  xdd = vector(mode="numeric",length =(2^l))
  for(id in 1:(2^l)){
    xdd[id] = fdis[id]*(Ydis[(p-1)*2^Lmax+id*2^(Lmax-l)+1]-Ydis[(p-1)*2^Lmax+(id-1)*2^(Lmax-l)+1]) - 0.5^(l+1)*(fdis[id]^2)
  }
  Gf[i]= sum(xdd)
  xdd = vector(mode="numeric",length =(2^(l-1)))
  for(id in 1:(2^(l-1))){
    xdd[id] = cdis[id]*(Ydis[(p-1)*2^Lmax+id*2^(Lmax-l+1)+1]-Ydis[(p-1)*2^Lmax+(id-1)*2^(Lmax-l+1)+1]) - 0.5^(l)*(cdis[id]^2)
  }
  Gc[i]= sum(xdd)
  }
  ## normalized weights function
  Wf = exp(Gf-max(Gf))
  Wc = exp(Gc-max(Gc))
  wf = Wf/sum(Wf)
  wc = Wc/sum(Wc)
  ## Effective sample size
  ESS1 = (sum(wf^2))^(-1)
  ## Compare with the threshold
  if(ESS1<=(N/2)){
  ## Maximal Resampling scheme
  minw = vector(mode = "numeric", length = N)
  for (i in 1:N) {
    minw[i] = min(wf[i],wc[i])
  }
  alpha = sum(minw)
  dice = runif(1, min = 0, max = 1)
  If = vector(mode = "numeric", length = N)
  Ic = vector(mode = "numeric", length = N)
  Zf = vector(mode = "numeric", length = N)
  Zc = vector(mode = "numeric", length = N)
  for (i in 1:N){
  Zf[i] = wf[i] - minw[i]
  Zc[i] = wc[i] - minw[i]
  }
  if ( dice <= alpha ){
    If <- rdiscrete(N, minw)
    Ic <- If
    } 
  if ( dice > alpha ){
    If <- rdiscrete(N, Zf)
    Ic <- rdiscrete(N, Zc)
    }
  for(i in 1:N){
    Xun1[p+1,i] = Xun[p+1,If[i]]
    Xunc1[p+1,i] = Xunc[p+1,Ic[i]]
  }
  }
  ## If Effective sample size larger than the threshold, we don't resample
  if(ESS1>(N/2)){
    Xun1[p+1,]=Xun[p+1,]
    Xunc1[p+1,]=Xunc[p+1,]
  }
}
## Estimates for Coupled Particle Filter, this is the summand term in the telescopic summation at l
for (i in 1:N) {
  Bb[i] = wf[i]*Xun[T+1,i] - wc[i]*Xunc[T+1,i]
}
Aa[l+1]=sum(Bb)
}
## Now, Aa include the summand for l from 1 to L , but l = 0 term has yet to be added into the map.

## Initialization
N = ceiling((2^L)*L) - 1
ko = vector(mode = "numeric",length = N)
Xun = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xun1 = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xun[1,]=ko
Xun1[1,]=ko
Wf = vector(mode = "numeric", length = N)
Gf = vector(mode = "numeric",length = N)
Bb = vector(mode = "numeric", length = N)
## Loop over time
for (p in 1:T) {
  ## Loop over particles 
  for(i in 1:N){
  ## Particle update for one time unit using Euler discretization
  ## Finer path particles
  fdis = vector(mode = "numeric", length =2)
  fdis[1] = Xun1[p,i]
  dw = rnorm(1, mean=0, sd = 1)
  fdis[2] = (-fdis[1]) + 0.5*dw + fdis[1]
  Xun[p+1,i]=fdis[2]
  ## G Function calcualtion
  xdd = vector(mode="numeric",length =(2^0))
  for(id in 1:(2^0)){
    xdd[id] = fdis[id]*(Ydis[(p-1)*2^Lmax+id*2^(Lmax)+1]-Ydis[(p-1)*2^Lmax+(id-1)*2^(Lmax)+1]) - 0.5^(1)*(fdis[id]^2)
  }
  Gf[i]= sum(xdd)
  }
  ## Normalized weights function
  Wf = exp(Gf-max(Gf))
  wf = Wf/sum(Wf)
  ## Effective Sample Size
  ESS2 = (sum(wf^2))^(-1)
  ## If ESS smaller or equal to the threshold, we resample
  if(ESS2<(N/2)){
  If = vector(mode = "numeric", length = N)
  If <- rdiscrete(N, Wf)
  for(i in 1:N){
    Xun1[p+1,i] = Xun[p+1,If[i]]
  }
  }
  ## If ESS larger than the threshold, we don't resample
  if(ESS2>=(N/2)){
    Xun1[p+1,]=Xun[p+1,]
  }
}
## Compute estimate value for l=0
for (i in 1:N) {
  Bb[i] = wf[i]*Xun[T+1,i]
}
Aa[1] = sum(Bb)
## Cauculate the final estimator value for MLPF
estMLPF[h] = sum(Aa)
}

## Cost&MSE for the plot
N = vector(length = L+1)
for(i in 0:L){
  N[i+1] = ceiling(2^(L-i)*L) - 1
}
K = vector(length = L+1)
K[1] = N[1]*T
for(i in 1:L){
  K[i+1] = N[i+1]*3*2^(i-1)*T
}
work_MLPF[L] = sum(K)
MSE_MLPF[L] = (mean(estMLPF)-tv)^2 + var(estMLPF)
var_MLPF[L] = var(estMLPF)
mean_MLPF[L] = mean(estMLPF)
}
date()
## nothing wrong with this, apparently
```


# OU MLPF Wasserstein
```{r}
# Gradient Estimator, dYt = Xt*dt+dBt, dXt = (-Xt)*dt + 0.5*dWt 
# Multilevel Particle Filter structure 
date()
mean1_PF = vector(mode = "numeric",length = 8)
MSE1_MLPF = vector(mode = "numeric",length = 8)
var1_MLPF = vector(mode = "numeric",length = 8)
work1_MLPF = vector(mode = "numeric",length = 8)
for(ltemp in 1:8){
T = 100
L = ltemp
R=200
estMLPF = vector(mode = "numeric",length = R)

for(h in 1:R){
Aa = vector(mode = "numeric", length = L+1)
for(l in 1:L){
N = 10*ceiling(2^(L-1.5*l))
ko = vector(mode = "numeric",length = N)
Xun = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xun1 = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xunc = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xunc1 = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
## replication of algorithm to get "true value"
Xun[1,]=ko
Xun1[1,]=ko
Xunc[1,]=ko
Xunc1[1,]=ko
wf = vector(mode = "numeric", length = N)
wc = vector(mode = "numeric", length = N)
Wf = vector(mode = "numeric", length = N)
Wc = vector(mode = "numeric", length = N)
Gf = vector(mode = "numeric",length = N)
Gc = vector(mode = "numeric",length = N)
Bb = vector(mode = "numeric", length = N)
for (p in 1:T) {
  for(i in 1:N){
  ## Particle update for one time unit
  dw = vector(mode = "numeric", length =(2^l))
  ## Finer path
  fdis = vector(mode = "numeric", length =(2^l+1))
  fdis[1] = Xun1[p,i]
  for (j in 1:(2^l)) {
      dw[j] = rnorm(1, mean=0, sd = sqrt(0.5^l))
      fdis[j+1] = (-fdis[j])*0.5^l + 0.5*dw[j] + fdis[j]
  }
  Xun[p+1,i]=fdis[(2^l+1)]
  ##Coarser path
  cdis = vector(mode = "numeric", length =(2^(l-1)+1))
  cdis[1] = Xunc1[p,i]
  for(j in 1:(2^(l-1))){
    cdis[j+1] = (-cdis[j])*0.5^(l-1) + 0.5*(dw[2*j]+dw[2*j-1]) + cdis[j]
  }
  Xunc[p+1,i]=cdis[(2^(l-1)+1)]
  ## G Function,note that there should only be one set of data,thus the y_{(k+1)*\Delta_{L}} - y_{(k)*\Delta_{L}} should be sampled only once
  xdd = vector(mode="numeric",length =(2^l))
  for(id in 1:(2^l)){
    xdd[id] = fdis[id]*(Ydis[(p-1)*2^Lmax+id*2^(Lmax-l)+1]-Ydis[(p-1)*2^Lmax+(id-1)*2^(Lmax-l)+1]) - 0.5^(l+1)*(fdis[id]^2)
  }
  Gf[i]= sum(xdd)
  xdd = vector(mode="numeric",length =(2^(l-1)))
  for(id in 1:(2^(l-1))){
    xdd[id] = cdis[id]*(Ydis[(p-1)*2^Lmax+id*2^(Lmax-l+1)+1]-Ydis[(p-1)*2^Lmax+(id-1)*2^(Lmax-l+1)+1]) - 0.5^(l)*(cdis[id]^2)
  }
  Gc[i]= sum(xdd)
  }
  ## normalized weights function
  Wf = exp(Gf-max(Gf))
  Wc = exp(Gc-max(Gc))
  wf = Wf/sum(Wf)
  wc = Wc/sum(Wc)
  ## Effective Sample Size
  ESS1 = (sum(wf^2))^(-1)
  ## If the ESS is smaller than the threshold, we resample
  if(ESS1<=(N/2)){
    If = vector(mode = "numeric",length = N)
    Ic = vector(mode = "numeric",length = N)
    ## Below is Wasserstein resampling method
    dice = runif(N,min = 0,max = 1)
    ## First we sort the particles
    uc1 = Xunc[p+1,]
    csort = uc1[order(uc1,decreasing = FALSE)]
    ## Sort weights the same way as the particles
    wcsort = wc[order(uc1,decreasing = FALSE)]
    ## Get the index
    cumc = cumsum(wcsort)
    Ic = findInterval(dice,cumc)+1
    ## Same step for finer path particles
    Uhat = Xun[p+1,]
    usort = Uhat[order(Uhat,decreasing = FALSE)]
    wsort = w[order(Uhat,decreasing = FALSE)]
    cum = cumsum(wsort)
    If = findInterval(dice,cum,rightmost.closed = TRUE)+1
    ## Update
    for(i in 1:N){
      Xun1[p+1,i] = usort[ix[i]]
      Xunc1[p+1,i] = csort[ix[i]]
    }
  }
  ## If the ESS is larger than the threshold, we don't resample
  if(ESS1>(N/2)){
    Xun1[p+1,]=Xun[p+1,]
    Xunc1[p+1,]=Xunc[p+1,]
  }
}
for (i in 1:N) {
  Bb[i] = wf[i]*Xun[T+1,i] - wc[i]*Xunc[T+1,i]
}
Aa[l+1]=sum(Bb)
}
## Now, Aa include the summand for l from 1 to L , but l = 0 term has yet to be added into the map.
N = ceiling((2^L)*L) - 1
ko = vector(mode = "numeric",length = N)
Xun = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xun1 = matrix(c(1:((T+1)*N)),nrow=(T+1),ncol=N)
Xun[1,]=ko
Xun1[1,]=ko
Wf = vector(mode = "numeric", length = N)
Gf = vector(mode = "numeric",length = N)
Bb = vector(mode = "numeric", length = N)
for (p in 1:T) {
  for(i in 1:N){
  ## Particle update for one time unit
  ## Finer path Euler update
  fdis = vector(mode = "numeric", length =2)
  fdis[1] = Xun1[p,i]
  dw = rnorm(1, mean=0, sd = 1)
  fdis[2] = (-fdis[1]) + 0.5*dw + fdis[1]
  Xun[p+1,i]=fdis[2]
  ## G Function calculation
  xdd = vector(mode="numeric",length =(2^0))
  for(id in 1:(2^0)){
    xdd[id] = fdis[id]*(Ydis[(p-1)*2^Lmax+id*2^(Lmax)+1]-Ydis[(p-1)*2^Lmax+(id-1)*2^(Lmax)+1]) - 0.5^(1)*(fdis[id]^2)
  }
  Gf[i]= sum(xdd)
  }
  ## Normalized weights function
  Wf = exp(Gf-max(Gf))
  wf = Wf/sum(Wf)
  ## ESS calculation
  ESS2 = (sum(wf^2))^(-1)
  ## If ESS smaller or equal to threshold, we resample
  if(ESS2<(N/2)){
    ## Getting index
    If = vector(mode = "numeric", length = N)
    If <- rdiscrete(N, Wf)
    ## Update
    for(i in 1:N){
      Xun1[p+1,i] = Xun[p+1,If[i]]
    }
  }
  ## If ESS smaller or equal to threshold, we resample
  if(ESS2>=(N/2)){
    Xun1[p+1,]=Xun[p+1,]
  }
}
for (i in 1:N) {
  Bb[i] = wf[i]*Xun[T+1,i]
}
Aa[1] = sum(Bb)
## Cauculate the final estimator value
estMLPF[h] = sum(Aa)
}

## Cost&MSE for the plot
N = vector(length = L+1)
for(i in 0:L){
  N[i+1] = ceiling(2^(L-i)*L) - 1
}
K = vector(length = L+1)
K[1] = N[1]*T
for(i in 1:L){
  K[i+1] = N[i+1]*3*2^(i-1)*T
}
work1_MLPF[L] = sum(K)
MSE1_MLPF[L] = (mean(estMLPF)-tv)^2 + var(estMLPF)
var1_MLPF[L] = var(estMLPF)
mean1_MLPF[L] = mean(estMLPF)
}
date()
## nothing wrong with this, apparently
```