[R] Is R the right choice for simulating first passage times of random walks?

[Paul Menzel]

Dear R folks,


Am Donnerstag, den 28.07.2011, 01:36 +0200 schrieb Paul Menzel:

> I need to simulate first passage times for iterated partial sums. The
> related papers are for example [1][2].
> 
> As a start I want to simulate how long a simple random walk stays
> negative, which should result that it behaves like n^(-½). My code looks
> like this.
> 
> -------- 8< -------- code -------- >8 --------
> n = 100000 # number of simulations
> 
> length = 100000 # length of iterated sum
> z = rep(0, times=length + 1)
> 
> for (i in 1:n) {
> 	x = c(0, sign(rnorm(length)))
> 	s = cumsum(x)
> 	for (i in 1:length) {
> 		if (s[i] < 0 && s[i + 1] >= 0) {
> 			z[i] = z[i] + 1
> 		}
> 	}
> }
> plot(1:length(z), z/n)
> curve(x**(-0.5), add = TRUE)
> -------- 8< -------- code -------- >8 --------

Of course the program above is not complete, because it only checks for
the first passage from negativ to positive. `if (s[2] < 0) {}` should be
added before the for loop.

> This code already runs for over half an hour on my system¹.
> 
> Reading about the for loop [3] it says to try to avoid loops and I
> probably should use a matrix where every row is a sample.
> 
> Now my first problem is that there is no matrix equivalent for
> `cumsum()`. Can I use matrices to avoid the for loop?

I mean the inner for loop.

Additionally I wonder if `cumsum` is really faster or if I should sum
the elements by myself and check after every step if the walk gets
non-negative/0. With a length of 1000000 this should save some cycles.
On the other hand adding numbers should be really fast and adding checks
in between could potentially be slower.

> My second question is, is R the right choice for such simulations? It
> would be great when R can also give me a confidence interval(?) and also
> try to fit a curve through the result and give me the rule of
> correspondence(?) [4]. Do you have pointers for those?
> 
> I glanced at simFrame [5] and read `?simulate` but could not understand
> it right away and think that this might be overkill.
> 
> Do you have any suggestions?


Thanks,

Paul


> ¹ AMD Athlon(tm) X2 Dual Core Processor BE-2350, 2,1 GHz
> 
> 
> [1] http://www-stat.stanford.edu/~amir/preprints/irw.ps
> [2] http://arxiv.org/abs/0911.5456
> [3] http://cran.r-project.org/doc/manuals/R-intro.html#Repetitive-execution
> [4] https://secure.wikimedia.org/wikipedia/en/wiki/Function_(mathematics)
> [5] http://finzi.psych.upenn.edu/R/library/simFrame/html/runSimulation.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 198 bytes
Desc: This is a digitally signed message part
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20110728/d8a1313d/attachment.bin>