Path of Ornstein Uhlenbeck Process

library(sde) # call sde for the simulation of SDE

## Loading required package: MASS
## Loading required package: stats4
## Loading required package: fda
## Loading required package: splines
## Loading required package: Matrix
## 
## Attaching package: 'fda'
## 
## The following object is masked from 'package:graphics':
## 
##     matplot
## 
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## sde 2.0.13
## Companion package to the book
## 'Simulation and Inference for Stochastic Differential Equations With R Examples'
## Iacus, Springer NY, (2008)
## To check the errata corrige of the book, type vignette("sde.errata")

rm(list = ls())     # clear objects  


# form the layaut for 3x2 graphics
layout(mat=matrix(1:6, 3, 2, byrow=FALSE))

# x values
# x=(l-N/2)/sqrt(N) where l is the number of filled sites
xax<-seq(-2,2,0.1) 

# Propagator obtained by path integral
propagator<-function(xi,t){
  pro<-sqrt(2/(pi*(1-exp(-4*t))))*exp(-2*(xi*exp(-2*t)-xax)^2/(1-exp(-4*t)))
  return(pro)
}


d <- expression(-2*x) # drift value to be used by sde.sim()
s <- expression(1)  # sigma value to be used by sde.sim()

x0 <- 0 # Initial x value, i.e., l=N/2
time <- 1000 # Horizon of simulation
numb <- 10000 # Number of simulation steps
par(cex.lab=1.5) # Put the axis font size

# X1, X2 and X3 keep the x values of the process
sde.sim(X0=x0,T=time,N=numb,drift=d,sigma=s) -> X1

## sigma.x not provided, attempting symbolic derivation.

plot(X1,xlab="t",ylab="path",cex.lab=2)
text(400, -1.5, bquote(x[0]==0), cex=1.5)

sde.sim(X0=5,T=time,N=numb,drift=d,sigma=s) -> X2

## sigma.x not provided, attempting symbolic derivation.

plot(X2,xlab="t",ylab="path")
text(20, 3, bquote(x[0]==5), cex=1.5)

sde.sim(X0=10,T=time,N=numb,drift=d,sigma=s) -> X3

## sigma.x not provided, attempting symbolic derivation.

plot(X3,xlab="t",ylab="path")
text(20, 6, bquote(x[0]==10), cex=1.5)

# Distribution of the process is plotted with histogram
# and compared with the analytic expression "propagator". 

hist(X1,xlim=c(-2,2),breaks=50,prob=TRUE,xlab="x",ylab="PDF",main="",cex.lab=1.5)
lines(seq(-2,2,0.1),propagator(x0,time),col="red",lwd=2)
text(-1, 0.6, bquote(analytic), cex=1.5,col="red")
text(0, 0.4, bquote(numeric), cex=1.5)

hist(X2,xlim=c(-2,2),breaks=50,prob=TRUE,xlab="x",ylab="PDF",main="",cex.lab=1.5)
lines(seq(-2,2,0.1),propagator(5,time),col="red",lwd=2)
text(0, 0.4, bquote(numeric), cex=1.5)
text(-1, 0.6, bquote(analytic), cex=1.5,col="red")

hist(X3,xlim=c(-2,2),breaks=50,prob=TRUE,xlab="x",ylab="PDF",main="",cex.lab=1.5)
lines(seq(-2,2,0.1),propagator(10,time),col="red",lwd=2)
text(0, 0.4, bquote(numeric), cex=1.5)
text(-1, 0.6, bquote(analytic), cex=1.5,col="red")

#dev.off()