Heston Simulation Monte Carlo: Slow R code

Question

I need to simulate the stock price, that follows stochastic volatility process (Heston Model). I already asked, how to speed up my loops, but for this case I´m not able to use some tips due to the V[i-1] dependence.

Basically the code is: V is the volatility of the stock and S is the stock price. And: a,b,c... constants.

Here is the code:

V[1] <- 0.04
S[1] <- 40
U <- matrix(NA, nrow=100000, ncol=200, byrow=TRUE)

### Function ###
Inv.Phi <- function(y){
              if (y <= p) {0} else {log(1-p)}
}

### Simulation ####
for(j in 1:100000){
  for(i in 2:200){
    m <- V[i-1] * c
    n <- V[i-1] * d
    phi <- n/m

    if(phi <= 1.5){ 
       Z <- rnorm(1)
       V[i] <- rnorm(1) * e
       K <- V[i-1] * f
    }else{ 
       p <- (phi-1) / (phi+1)
       u <- runif(1)
       V[i] <- Inv.Phi(u)
       K <- V[i-1] * g
    }
 S[i] <- S[i-1] * exp(K * V[i-1]) * exp(V[i] * rnorm(1))
 }
U[j,] = S
}

Any suggestion to speed up this process! I known, I´m using a lot of bad things for R, but I couldn´t figure out a better solution.

How exactly is this different than your other question? http://stackoverflow.com/questions/15534270/stock-price-simulation-r-code-slow-monte-carlo — Paul Hiemstra, Mar 22 '13 at 21:03
Is anyone else getting the suspicion we are asked to do this person's business school homework? — IRTFM, Mar 22 '13 at 22:20
If it were a [Heston model](http://en.wikipedia.org/wiki/Heston_model) (it looks different: there is a threshold in your model), you could try `sde::sde.sim` to simulate the volatility, and vectorize your code, as in your previous question, to compute the prices. If this is not sufficient, you can try to rewrite the loop in C/C++, e.g., with [Rcpp](https://github.com/hadley/devtools/wiki/Rcpp). — Vincent Zoonekynd, Mar 22 '13 at 22:34
As I see I wasn´t supposed to open another question. I apologize for that. But I´m not able to put inside the matrix the if else and the initial price. And also due to the fact that the stock price depends also to another stochastic factor. Anyway sorry for the inconvenient — PereMkB, Mar 22 '13 at 22:53
I'm sure it would be an interesting exercise to *speed up/fix* your homework, bu as it stands we can't. You haven't defined `c` or `d` or `e` or `f` and indeed `g`! — Simon O'Hanlon, Mar 22 '13 at 22:54
Hi, all the constants are define before. This is just part of the code. I´ve just tried to give you a picture of the problem. These constants are huge, and depend in a lot of inputs. — PereMkB, Mar 22 '13 at 22:58
And again, just to make sure, I really appreciate all help from my other topic. I´m still trying to implement in this case and with jump process also. I´m truly sorry if pass the wrong idea. @Vicent, I used the Andersen Scheme (not Euler Scheme) to discretize the stochastic process. Thanks for the tip — PereMkB, Mar 22 '13 at 23:02

Ben Bolker · Accepted Answer · 2013-03-23T14:17:51.657

0

You will probably speed things up a lot if you get rid of the outer loop and construct an entire column of your matrix at once:

nsim <- 1e5
nt <- 200
U <- matrix(NA, nrow=nsim, ncol=200, byrow=TRUE)
V_last <- V <- rep(0.04, nsim)
U[,1] <- 40

Inv.phi <- function(y,p) ifelse (y <= p, 0, log(1-p))
for(i in 2:nt) {
    m <- V_last * c
    n <- V_last * d
    phi <- n/m ## ??? as currently stated this is just a constant d/c ??
               ## presumably there is a typo somewhere??
    bt <- (phi<=1.5) ## below-threshold
    V[bt] <- rnorm(length(bt)) * e
    V[!bt] <- Inv.phi(runif(length(!bt)), (phi[!bt]-1) / (phi[!bt]+1))
    K <- V_last*ifelse(bt,f,g)
    U[,i] <- U[,i-1] * exp(K * V_last) * exp(V * rnorm(nsim))
    V_last <- V
}

I think this works but can't test it because you haven't given a reproducible example ...

edited Mar 23 '13 at 14:17

answered Mar 23 '13 at 00:49

Ben Bolker

211,554
25
370
453

thanks a lot for your attention! When I was shrinking the code, I forgot to delete that “Z”. Yes, d and c are constants. I´ve just tried your idea and works much, much faster than with both loops. However, there is a warning regarding the if/else. It warns that the size of phi>1, so it is using only the first one. When I´m running the complete simulation, the expected answer (for an example) is approximately 4.47, but now is resulting near 4,58. I´ll try to find out what is happening. Thanks a lot again – PereMkB Mar 23 '13 at 09:22
fixed to allow for vector `phi`. The issue is not that `c` and `d` are constants, it's that `phi = (V_last*d)/(V_last*c) = d/c`, so `phi` will also always be a constant ... – Ben Bolker Mar 23 '13 at 14:19

Heston Simulation Monte Carlo: Slow R code

1 Answers1