Search code examples
rgraphsimulation

Brownian motion simulation in R

I have this process š‘‹š‘”=-3š‘”+2šµš‘” that I want to simulate using R. šµš‘” is a standard brownian motion.

However, I have figured that š‘‹š‘” is not a brownian motion, since its mean is š”¼[š‘‹š‘”]=š”¼[-3š‘”+2šµš‘”]=-3š‘”+š”¼[2šµš‘”]=-3š‘” (not 0) and the variance is š”¼[(2šµš‘”)^2]=4š”¼[šµ^2š‘”]=4š‘”.

When looking online on how to simulate brownian motion, I see that usually normal distribution is used, for example rnorm(n,0,t) my question is, how would I simulate this š‘‹š‘” when it is not a standard brownian motion?

Solution

  • It is perfectly valid to combine a linear term (or any other function) with a random walk term. Here is a basic skeleton...

    x <- 1:20
xt <- function (x) { .5*x +                       # linear term plus 
                      cumsum(rnorm(length(x)))    # random walk term
                   }   
                                                    
#xt(x) # display one instance if desired
matplot(rbind(0,replicate(100,xt(x))), type='l') # start at origin w/ rbind of 0's

    You, of course, want xt to look something like:

    xt <- function (x) { (-3)*x+cumsum(rnorm(length(x)))*2} # linear term plus 
                                                        # random walk term

    Re. your 2nd question (below) You mean 2, not 4 (that's the variance value, not sd). It's the same either way (except for floating point rounding error hence the use of all.equal() here).

    xt <- function (x) { (-3)*x+cumsum(rnorm(length(x))*2)}
xt2 <- function (x) { (-3)*x+cumsum(rnorm(length(x),sd=2))}
set.seed(123L)
a <- xt(x)
set.seed(123L)
b <- xt2(x)
all.equal(a,b)

    Added: One last thing I can mention is matplot() is a WONDERFUL tool for displaying the results of Monte Carlo simulations like random walks.

    Added #2: Re. Brownian motion rather than simple gaussian random walks, here is a quickie script.

    x <- 1000    
xbt_gauss <- function (x) {cumsum(rnorm(x))}
# xbt_gauss(x) # examine 1 instance

plot(NULL,
     xlim=c(-100,100),
     ylim=c(-100,100))

for (i in 1:50) lines(cbind(xbt_gauss(x),xbt_gauss(x)),
                      type='l',
                      col=i)

    Note: this is a reasonable use of a for loop in R (one of few, btw). One could also use foreach() in parallel coding if one wanted to do massive models here.

    Of course one could include additive or other terms here as your model dictates.