4 ms·
your point at the end is good - if we're only interested in the final price, then we should be able to do this in constant memory, independent of the number of
by shoo 13y ago
your point at the end is good - if we're only interested in the final price, then we should be able to do this in constant memory, independent of the number of timesteps and samples.
edit: on my machine the running time seems dominated by the generation of random numbers, rather than anything else.
perhaps another optimisation is to elimate the `exp` from the inner loop, and just do addition. if you're only interested in the final price, you can do one `exp` at the end.
i.e. change
paths[1, i] = st = s0
for j in 2:M
st *= exp(a + b*randn())
paths[j, i] = st
end
to this
paths[1, i] = log_st = log(s0)
for j in 2:M
log_st += (a + b*randn())
paths[j, i] = log_st
end
- Retric 13y agoOne solution to this is use a circular data structure filled with random numbers. Then dedicate a thread to updating that ring with new random numbers. Meanwhile you pull random numbers as fast as you want. It's not good for encryption, but for simulation it's fine. http://en.wikipedia.org/wiki/Circular_buffer http://en.wikipedia.org/wiki/Circular_buffer Another approach is to keep reseeding a fast RNG with a slow but higher quality RNG.
- StefanKarpinski 13y agoYou don't have to do anything quite that complicated. You can just pre-generate an array of Gaussian random samples and do table lookup in the inner loop: paths = zeros(Float64, M, I) rands = zeros(Float64, M-1, 1) for i in 1:I randn!(rands) paths[1, i] = st = s0 for j in 2:M st *= exp(a + b*rands[j-1]) paths[j, i] = st end end That shaves off another 35% percent of runtime on my machine.
- jwmerrill 13y agoNice! In case it isn't obvious, randn! takes an existing vector and fills it up with new random normal variates. This allows computing a bunch of random variates at once, which is a little more efficient, without having to allocate a new vector every time through the loop. Another option would be to generate all M-1xI random numbers ahead of time, but that would be less memory efficient.
- StefanKarpinski 13y agoI tried that too, but it wasn't any faster – you've already gotten most of the gains to be had by lifting the random number generation out of the innermost loop. So I went with the option that uses less memory.
- shoo 13y agoi tried the approach of generating them ahead of time and then loading them from disk. it only uses 80mb of memory for this example. the difference i see between generating numbers on-the-fly using boost and reading them from disk is 6.41s vs 0.32s, ie a 20x speedup, for 10 repetitions of 100,000 paths over 100 timesteps. edit: here's the code - https://gist.github.com/fcostin/3ed3512ad167adf9c785 https://gist.github.com/fcostin/3ed3512ad167adf9c785
- jwmerrill 13y agoI would not recommend that circular buffer idea for quantitative calculations. Consider averaging n numbers from a buffer that you've filled with m random variates with m < n: the standard deviation of your average will decrease like 1/sqrt(m), not 1/sqrt(n), so you might as well have just averaged the independent random numbers you produced instead of reusing some of them. Maybe the circular buffer idea would be useful in non-quantitative contexts, though, e.g. for simulating something in a game or a movie, where all you need is the appearance of randomness.
- StefanKarpinski 13y agoAgreed, that circular buffer thing is bound to do one of two things based on whether the RNG thread or the work thread is faster: a) work slower than RNG – do too much work, wasting CPU; b) work faster than RNG – produce very poor randomness. Neither of these situations is good, and dealing with two threads in the first place is a bit of a complication. If you're going to do that, it's better to synchronize the two threads so that the work thread blocks when it needs more randomness than the RNG thread has produced and vice versa.