Reputation: 13
I have a problem with a nested loop, here's the code:
T=2
dt=0.001
K=floor(T/dt)
t= seq(0,T,dt)
BM<- matrix(0,2000,1000)
for (i in 1:1000){
for (k in 1:K-1){
BM[k+1,i]=BM[k,i]+sqrt(dt)*rnorm(1)
}
}
the error message is:
Error in BM[k + 1, i] = BM[k, i] + sqrt(dt) * rnorm(1) :
replacement has length zero
What is wrong in the loop? The goal of it is to simulate 1000 sample paths of Brownian motion.
Upvotes: 1
Views: 2377
Reputation: 61154
Here's a fully vectorized approach which gives the same result and it is faster than using for
loops
B <- mat.or.vec(K,1000) # pre-allocate
set.seed(1)
B[-1, ] <- matrix(rnorm((K-1)*100),K-1,1000) * sqrt(dt)
B <- apply(B, 2, cumsum)
Vectorized functions speed up your code as you can see in the benchmark results below:
test elapsed relative
1 vectorized() 2.070 1.000
2 for.loop() 181.768 87.811
The vectorized approach is about 88 times faster than your two nested for
-loops. Codes for the benchmark are shown below.
for.loop <- function(){ # this is your approach
BM<- matrix(0,2000,1000)
set.seed(1)
for (i in 1:1000){
for (k in 1:(K-1)){
BM[k+1, i]=BM[k,i]+sqrt(dt)*rnorm(1)
}
}
return(BM)
}
vectorized <- function(){ # this is my approach
B <- mat.or.vec(K,1000) # pre-allocate
set.seed(1)
B[-1, ] <- matrix(rnorm((K-1)*100),K-1,1000) * sqrt(dt)
B <- apply(B, 2, cumsum)
}
library("rbenchmark")
benchmark(vectorized(),
for.loop(),
replications=10,
columns=c('test', 'elapsed', 'relative'),
order = "relative")
Comparing the first 5 rows and 5 cols to check results:
> vectorized()[1:5, 1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00000000 0.000000000 0.000000000 0.000000000 0.000000000
[2,] -0.01981021 -0.009865464 -0.008371525 0.015834886 0.001562361
[3,] -0.01400290 -0.037887974 -0.037253019 0.010349073 -0.024924006
[4,] -0.04042779 -0.098675011 -0.073133176 0.018483702 -0.066867547
[5,] 0.01001941 -0.047455576 -0.048955758 -0.001085329 -0.123356698
> for.loop()[1:5, 1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00000000 0.000000000 0.000000000 0.000000000 0.000000000
[2,] -0.01981021 -0.009865464 -0.008371525 0.015834886 0.001562361
[3,] -0.01400290 -0.037887974 -0.037253019 0.010349073 -0.024924006
[4,] -0.04042779 -0.098675011 -0.073133176 0.018483702 -0.066867547
[5,] 0.01001941 -0.047455576 -0.048955758 -0.001085329 -0.123356698
Upvotes: 3
Reputation: 42629
This is the problem:
for (k in 1:K-1){
Your precedence is off. This will work:
for (k in 1:(K-1)){
Equivalently (for K > 0
):
for (k in seq(K-1)){
Vectorization is better yet:
for (k in seq(K-1)) {
BM[k+1,] <- BM[k,]+sqrt(dt)*rnorm(2000)
}
Upvotes: 2