Brownian Motion Trajectories

rm(list=ls())
BMtrend<-function(T=1,N=N,start.p=start.p,mu=mu,sigma=sigma){
  y.seq<-array(NA,N)
  y.seq[1]=start.p
  for(i in 2:N){
    y.seq[i]= y.seq[i-1]+  mu + sigma*rnorm(1)
  }
  y.seq
}

BM without trend

mu<-0 # here is no trend
sigma<-1
N<-1000
start.p<-0
y.seq.array<-NULL
for(i in 1:N){
  y.seq<-BMtrend(T=1,N=N,start.p=start.p,mu=mu,sigma=sigma)
  y.seq.array<-rbind(y.seq.array,y.seq)
}

par(mfrow=c(1,2))
plot(y.seq.array[1,],type="l",ylim=range(y.seq.array), xlab = "time", ylab = "phenotype",main= expression(paste("BM without Trend: ",paste(y[t]," = ",sigma," ",W[t]) )))
for(i in 2:100){
  points(1:N,y.seq.array[i,],type="l")
}

BM with trend

mu<-0.2 # here is the trend
sigma<-1
start.p<-0
y.seq.array<-NULL
for(i in 1:N){
  y.seq<-BMtrend(T=1,N=N,start.p=start.p,mu=mu,sigma=sigma)
  y.seq.array<-rbind(y.seq.array,y.seq)
}
plot(y.seq.array[1,],type="l",ylim=range(y.seq.array), xlab = "time", ylab = "phenotype",main= expression(paste("BM with Trend: ",paste(y[t]," = ",mu,t," + ",sigma,W[t]) )))
for(i in 2:N){
  points(1:N,y.seq.array[i,],type="l")
}