Tristan Linke
2011-Jul-02 23:07 UTC
[R] Simulating inhomogeneous Poisson process without loop
Dear all I want to simulate a stochastic jump variance process where N is Bernoulli with intensity lambda0 + lambda1*Vt. lambda0 is constant and lambda1 can be interpreted as a regression coefficient on the current variance level Vt. J is a scaling factor How can I rewrite this avoiding the loop structure which is very time-consuming for long simulations? for (i in 1:N){ ... N <- rbinom(n=1, size=1, prob=(lambda0+lambda1*Vt)) Vt <- ... + J*N .. } P.S. This is going towards the Duffie, Pan, Singleton 2000 Transform Pricing paper, here stochastic volatility with state-dependent correlated jumps (Eraker 2004). Thanks a lot in advance. [[alternative HTML version deleted]]
Ben Bolker
2011-Jul-05 11:53 UTC
[R] Simulating inhomogeneous Poisson process without loop
Tristan Linke <tristan.linke <at> gmail.com> writes:> > Dear all > > I want to simulate a stochastic jump variance process where N is Bernoulli > with intensity lambda0 + lambda1*Vt. lambda0 is constant and lambda1 can be > interpreted as a regression coefficient on the current variance level Vt. J > is a scaling factor > > How can I rewrite this avoiding the loop structure which is very > time-consuming for long simulations? > > for (i in 1:N){ > ... > N <- rbinom(n=1, size=1, prob=(lambda0+lambda1*Vt)) > Vt <- ... + J*N > .. > }Is it a typo that you are using N both as your loop variable and as part of the state of your simulation?> P.S. This is going towards the Duffie, Pan, Singleton 2000 Transform Pricing > paper, here stochastic volatility with state-dependent correlated jumps > (Eraker 2004).I don't think there's any way to rearrange this completely without loops, because each call of rbinom() depends on the previous updating. One thing that might speed things up a lot would be to pick N random uniform variates as a single vector in advance: rvec <- runif(N) for (i in 1:N) { Vt <- ... + J*(rvec[i]<lambda0+lambda1*Vt) } Depending on how badly you need it to be fast/how much effort you're willing to invest, you could try compiling/using inline or Rcpp etc. The other point is that if the probability stays constant over several steps when the jump doesn't happen, and if the probability is low, you could skip a few loop steps by drawing a geometric RV to determine the time of the next jump, rather than testing one step at a time.