bm ou eb path

rm(list=ls())
bm1d<-function(params,x=x,sims=sims){
  sigma<-params["sigma"]
  for(i in 2:sims){
    x[i]<-x[i-1] + sigma*rnorm(1)
  }
  return(x)
}

set.seed(668768)

sims<-266

x1<-array(NA,sims)
sigma<-0.1
x1[1]<-0
params<-c(sigma)
names(params)<-c("sigma")

x11<-bm1d(params,x=x1,sims=sims)
plot(x11,type="l",xlab="",ylab="",main=bquote(sigma~"="~ 0.1))
points(1,0,pch=5,col="blue")

#x2
#0.25
x2<-array(NA,sims)
sigma<-0.25
x2[1]<-0
params<-c(sigma)
names(params)<-c("sigma")

x22<-bm1d(params,x=x2,sims=sims)
plot(x22,type="l",xlab="",ylab="",main=bquote(sigma~"="~ 0.25))
points(1,0,pch=5,col="blue")

#x3
#0.8
x3<-array(NA,sims)
sigma<-0.8
x3[1]<-0
params<-c(sigma)
names(params)<-c("sigma")

x33<-bm1d(params,x=x3,sims=sims)
plot(x33,type="l",xlab="",ylab="",main=bquote(sigma~"="~ 0.8))
points(1,0,pch=5,col="blue")

#x4
#1.2
x4<-array(NA,sims)
sigma<-0.8
x4[1]<-0
params<-c(sigma)
names(params)<-c("sigma")

x44<-bm1d(params,x=x4,sims=sims)
plot(x44,type="l",xlab="",ylab="",main=bquote(sigma~"="~ 1.2))
points(1,0,pch=5,col="blue")

## OU


ou1d<-function(params,x=x,sims=sims){
  alpha<-params["alpha"]
  theta<-params["theta"]
  sigma<-params["sigma"]
  for(i in 2:sims){
    x[i]<-x[i-1] +  alpha*(theta-x[i-1])+ sigma*rnorm(1)}
  return(x)
}


set.seed(668768)

sims<-266

x<-array(NA,sims)
alpha<-1e-3
theta<-22
sigma<-1
x[1]<-0
params<-c(alpha,theta,sigma)
names(params)<-c("alpha","theta","sigma")


x1<-ou1d(params,x=x,sims=sims)
plot(x1,type="l",xlab="",ylab="",main=bquote(alpha~"="~ 0.00001~","~theta~"="~22))
points(1,0,pch=5,col="blue")
abline(h=theta,col="red")

#

y<-array(NA,sims)
alpha<-0.12
theta<-22
sigma<-1
y[1]<-0
params<-c(alpha,theta,sigma)
names(params)<-c("alpha","theta","sigma")

y1<-ou1d(params,x=y,sims=sims)
plot(y1,type="l",xlab="",ylab="",main=bquote(alpha~"="~ 0.8~","~theta~"="~22))
points(1,0,pch=5,col="blue")
abline(h=theta,col="red")

#

x<-array(NA,sims)
alpha<-1e-3
theta<--22
sigma<-1
x[1]<-0
params<-c(alpha,theta,sigma)
names(params)<-c("alpha","theta","sigma")

x2<-ou1d(params,x=x,sims=sims)

plot(x2,type="l",xlab="",ylab="",main=bquote(alpha~"="~ 0.00001~","~theta~"="~-22))
points(1,0,pch=5,col="blue")

#

y<-array(NA,sims)
alpha<-0.12
theta<- -22
sigma<-1
y[1]<-0
params<-c(alpha,theta,sigma)
names(params)<-c("alpha","theta","sigma")


y2<-ou1d(params,x=y,sims=sims)
plot(y2,type="l",xlab="",ylab="",main=bquote(alpha~"="~ 0.8~","~theta~"="~-22))
points(1,0,pch=5,col="blue")
abline(h=theta,col="red",lty=2,lwd=2)

## EB


eb1d<-function(params,x=x,sims=sims){
  r<-params["r"]
  sigma<-params["sigma"]
  for(i in 2:sims){
    x[i]<-x[i-1] +  exp(r*i/2)*sigma*rnorm(1)
  }
  return(x)
}


set.seed(668768)

sims<-266

xx1<-array(NA,sims)
r<--0.02
sigma<-1
xx1[1]<-0
params<-c(r,sigma)
names(params)<-c("r","sigma")


xx1<-eb1d(params,x=xx1,sims=sims)
plot(xx1,type="l",xlab="",ylab="",main=bquote(r~"="~ -0.02~","~sigma~"="~1))
points(1,0,pch=5,col="blue")

#
yy1<-array(NA,sims)
r<--0.025
sigma<-1
yy1[1]<-0
params<-c(r,sigma)
names(params)<-c("r","sigma")


yy1<-eb1d(params,x=yy1,sims=sims)
plot(yy1,type="l",xlab="",ylab="",main=bquote(r~"="~ -0.025~","~sigma~"="~1))
points(1,0,pch=5,col="blue")

#

xx2<-array(NA,sims)
r<--0.03
sigma<-1
xx2[1]<-0
params<-c(r,sigma)
names(params)<-c("r","sigma")

xx2<-eb1d(params,x=xx2,sims=sims)
plot(xx2,type="l",xlab="",ylab="",main=bquote(r~"="~ -0.03~","~sigma~"="~1))
points(1,0,pch=5,col="blue")

#

yy2<-array(NA,sims)
r<--0.035
sigma<-1
yy2[1]<-0
params<-c(r,sigma)
names(params)<-c("r","sigma")


yy2<-eb1d(params,x=yy2,sims=sims)
plot(yy2,type="l",xlab="",ylab="",main=bquote(r~"="~ -0.035~","~sigma~"="~1))
points(1,0,pch=5,col="blue")

###########################################
###########Elaborate Figures###############
###########################################


# https://stackoverflow.com/questions/66077931/r-how-to-use-proper-minus-symbol-for-negative-axes-in-plot


# plot(-3:3, -3:3, xaxt="n", yaxt="n", xlab = "\u22123:3", ylab = "\u22123:3")
# 
# xat <- axTicks(1, usr=par("usr")[1:2])
# labs <- gsub("-", "\U2212", print.default(xat))
# axis(1, at=xat, labels=labs)
# yat <- axTicks(2, usr=par("usr")[1:2])
# labs <- gsub("-", "\U2212", print.default(yat))
# axis(2, at=xat, labels=labs)


par(mfrow=c(1,3))


# bm
par(mar=c(4,4,2,2))
ylim<-range(x11,x22,x33,x44)

plot(NA,type="l",xlab="Time",ylab="Trait values",xaxt="n", yaxt="n",ylim=ylim, xlim=c(0,length(x11)),lwd = 2,main="Brownian Motion")

#plot(NA,type="l",xlab="Time",ylab="Trait values",ylim=ylim, xlim=c(0,length(x11)),lwd = 2,main="Brownian Motion")
points(1,0,pch=5,col="blue")

xat <- axTicks(1, usr=par("usr")[1:2])
labs <- gsub("-", "\U2212", print.default(xat))

## [1]   0  50 100 150 200 250

axis(1, at=xat, labels=labs)
yat <- axTicks(2, usr=par("usr")[1:2])

labs <- gsub("-", "\U2212", print.default(yat))

## [1] -20 -15 -10  -5   0   5

axis(2, at=yat, labels=labs)

lines(x11, lty=1 , lwd = 2 , col ="black")

lines(x22, lty=1 , lwd = 2 , col ="dodgerblue")

lines(x33, lty=1 , lwd = 2 , col ="green4")

lines(x44, lty=1 , lwd = 2 , col ="orange")

text.legend=c(expression(paste(sigma~"="~ 0.10))
              ,expression(paste(sigma~"="~ 0.25))
              ,expression(paste(sigma~"="~ 0.80))     
              ,expression( paste(sigma~"="~ 1.20)))

legend(x = 110, y = -16,legend=text.legend , lty=c(1 ,1),cex=0.75 , lwd=c(2 ,2) , col=c(" black" ,"dodgerblue","green4","orange"))

# ou

par(mar=c(4,4,2,2))
ylim<-range(x1,x2,y1,y2)

#plot(x1,type="l",xlab="Time",ylab="Trait values",ylim=ylim, lwd = 2,main="Ornstein--Uhlenbeck Process")

plot(NA,type="l",xlab="Time",ylab="Trait values",xaxt="n", yaxt="n",ylim=ylim, xlim=c(0,length(x11)),lwd = 2,main="Ornstein--Uhlenbeck Process")

points(1,0,pch=5,col="blue")

xat <- axTicks(1, usr=par("usr")[1:2])
labs <- gsub("-", "\U2212", print.default(xat))

## [1]   0  50 100 150 200 250

axis(1, at=xat, labels=labs)
yat <- axTicks(2, usr=par("usr")[1:2])

labs <- gsub("-", "\U2212", print.default(yat))

## [1] -30 -20 -10   0  10  20  30

axis(2, at=yat, labels=labs)

abline(h=c(-22,22),col="red",lty=2,lwd=2)

lines(x1, lty=1 , lwd = 2 , col ="black")

lines(y1, lty=1 , lwd = 2 , col ="dodgerblue")

lines(x2, lty=1 , lwd = 2 , col ="green4")

lines(y2, lty=1 , lwd = 2 , col ="orange")


text.legend=c(expression(paste(alpha~"="~ 0.120~","~theta~"="~-22~","~sigma~"="~1))
              ,expression(paste(alpha~"="~ 0.001~","~theta~"="~-22~","~sigma~"="~1))
              ,expression(paste(alpha~"="~ 0.120~","~theta~"="~22~","~sigma~"="~1))  
              ,expression( paste(alpha~"="~ 0.001~","~theta~"="~22~","~sigma~"="~1)))


legend(x = 100, y = 1 ,legend=text.legend , lty=c(1 ,1) ,cex=0.75, lwd=c(2 ,2) , col=c(" black" ,"dodgerblue","green4","orange"))



# eb

par(mar=c(4,4,2,2))

ylim<-range(xx1,xx2,yy1,yy2)

plot(NA,type="l",xlab="Time",ylab="Trait values",xaxt="n", yaxt="n",ylim=ylim, xlim=c(0,length(x11)),lwd = 2,main="Early Burst Process")

points(1,0,pch=5,col="blue")
# abline(h=c(-22,22),col="red",lty=2,lwd=2)

xat <- axTicks(1, usr=par("usr")[1:2])
labs <- gsub("-", "\U2212", print.default(xat))

## [1]   0  50 100 150 200 250

axis(1, at=xat, labels=labs)
yat <- axTicks(2, usr=par("usr")[1:2])

labs <- gsub("-", "\U2212", print.default(yat))

## [1] -10  -8  -6  -4  -2   0

axis(2, at=yat, labels=labs)


lines(xx1, lty=1 , lwd = 2 , col ="black")

lines(yy1, lty=1 , lwd = 2 , col ="dodgerblue")

lines(xx2, lty=1 , lwd = 2 , col ="green4")

lines(yy2, lty=1 , lwd = 2 , col ="orange")


text.legend=c(expression(paste(r~"="~ -0.020~","~sigma~"="~1))
              ,expression(paste(r~"="~ -0.025~","~sigma~"="~1))
              ,expression(paste(r~"="~ -0.030~","~sigma~"="~1))     
              ,expression(paste(r~"="~ -0.035~","~sigma~"="~1))     
)

legend(x = 150, y = -4 ,legend=text.legend , lty=c(1 ,1) ,cex=0.75, lwd=c(2 ,2) , col=c(" black" ,"dodgerblue","green4","orange"))