3

I need to perform a stock price simulation using R code. The problem is that the code is a little bit slow. Basically I need to simulate the stock price for each time step (daily) and store it in a matrix.

An example assuming the stock process is Geometric Brownian Motion

for(j in 1:100000){
    for(i in 1:252){
        S[i] <- S[i-1]*exp((r-v^2/2)*dt+v*sqrt(dt)*rnorm(1))
    }
    U[j,] <- S
}

Any suggestion to improve and speed up the code?

Dason
  • 60,663
  • 9
  • 131
  • 148
PereMkB
  • 29
  • 2
  • 5
  • 5
    Please provide a [reproducible example](http://stackoverflow.com/q/5963269). Your second for loop will throw an error because `S[i-1]` has zero-length when `i=1`. – Joshua Ulrich Mar 20 '13 at 20:58
  • Rewrite your code with `Rcpp` and you will get a speed-up of a few orders of magnitude. – Gary Weissman Mar 21 '13 at 01:52
  • Before going to `Rcpp` I would first try to get rid of the double for loop, and go for vectorisation. – Paul Hiemstra Mar 21 '13 at 06:35
  • 1
    In addition, you are growing the S and U vector inside the loop, this is very slow. Either preallocate these vectors, or use a solution like @Ferdinand.kraft suggested. – Paul Hiemstra Mar 21 '13 at 06:41

1 Answers1

2

Assuming S[0] = 1, you can build U as a follows:

Ncols <- 252

Nrows <- 100000

U <- matrix(exp((r-v^2/2)*dt+v*sqrt(dt)*rnorm(Ncols*Nrows)), ncol=Ncols, nrow=Nrows)

U <- do.call(rbind, lapply(1:Nrows, function(j)cumprod(U[j,])))

EDIT: using Joshua's and Ben's suggestions:

product version:

U <- matrix(exp((r-v^2/2)*dt+v*sqrt(dt)*rnorm(Ncols*Nrows)), ncol=Ncols, nrow=Nrows)

U <- t(apply(U, 1, cumprod))

sum version:

V <- matrix((r-v^2/2)*dt+v*sqrt(dt)*rnorm(Ncols*Nrows), ncol=Ncols, nrow=Nrows)

V <- exp( t(apply(V, 1, cumsum)) )

EDIT: as suggested by @Paul:

Execution time for each proposal (using 10000 rows instead of 10^5):

Using apply + cumprod

 user  system elapsed 
0.61    0.01    0.62

Using apply + cumsum

 user  system elapsed 
0.61    0.02    0.63

Using OP's original code

 user  system elapsed 
67.38    0.00   67.52

Notes: The times shown above are the third measures of system.time. The first two measures for each code were discarded. I've used r <- sqrt(2), v <- sqrt(3) and dt <- pi. In his original code, I've also replaced S[i-1] for ifelse(i==1,1,S[i-1]), and preallocated U.

Ferdinand.kraft
  • 12,579
  • 10
  • 47
  • 69