Код R для моделирования цены акций - Медленный - Монте-Карло

Мне нужно выполнить моделирование курса акций с использованием кода R. Проблема в том, что код немного медленный. В основном мне нужно моделировать курс акций для каждого временного шага (ежедневно) и сохранять его в матрице.

Пример, предполагающий, что складской процесс - это геометрическое броуновское движение.

for(j in 1:100000){
    for(i in 1:252){
        S[i] <- S[i-1]*exp((r-v^2/2)*dt+v*sqrt(dt)*rnorm(1))
    }
    U[j,] <- S
}

Есть предложения по улучшению и ускорению кода?


r montecarlo
person PereMkB    schedule 20.03.2013    source источник
comment
Приведите воспроизводимый пример. Второй цикл for выдаст ошибку, потому что S[i-1] имеет нулевую длину, когда i=1.   -  person Joshua Ulrich    schedule 21.03.2013
comment
Перепишите свой код с Rcpp, и вы получите ускорение на несколько порядков.   -  person Gary Weissman    schedule 21.03.2013
comment
Прежде чем перейти к Rcpp, я бы сначала попытался избавиться от двойного цикла for и перейти к векторизации.   -  person Paul Hiemstra    schedule 21.03.2013
comment
Кроме того, вы увеличиваете вектор S и U внутри цикла, это очень медленно. Либо предварительно распределите эти векторы, либо используйте решение вроде @ Ferdinand.kraft.   -  person Paul Hiemstra    schedule 21.03.2013

Ответы (1)


arrow_upward
2
arrow_downward

Предполагая S[0] = 1, вы можете построить U следующим образом:

Ncols <- 252

Nrows <- 100000

U <- matrix(exp((r-v^2/2)*dt+v*sqrt(dt)*rnorm(Ncols*Nrows)), ncol=Ncols, nrow=Nrows)

U <- do.call(rbind, lapply(1:Nrows, function(j)cumprod(U[j,])))

РЕДАКТИРОВАТЬ: используя предложения Джошуа и Бена:

версия продукта:

U <- matrix(exp((r-v^2/2)*dt+v*sqrt(dt)*rnorm(Ncols*Nrows)), ncol=Ncols, nrow=Nrows)

U <- t(apply(U, 1, cumprod))

Суммарная версия:

V <- matrix((r-v^2/2)*dt+v*sqrt(dt)*rnorm(Ncols*Nrows), ncol=Ncols, nrow=Nrows)

V <- exp( t(apply(V, 1, cumsum)) )

РЕДАКТИРОВАТЬ: как было предложено @Paul:

Время выполнения каждого предложения (с использованием 10000 строк вместо 10 ^ 5):

Использование apply + cumprod

 user  system elapsed 
0.61    0.01    0.62

Использование apply + cumsum

 user  system elapsed 
0.61    0.02    0.63

Использование исходного кода OP

 user  system elapsed 
67.38    0.00   67.52

Примечания: Время, показанное выше, является третьим тактом system.time. Первые две меры для каждого кода были отброшены. Я использовал r <- sqrt(2), v <- sqrt(3) и dt <- pi. В его исходном коде я также заменил S[i-1] на ifelse(i==1,1,S[i-1]) и предварительно выделил U.

person Ferdinand.kraft    schedule 20.03.2013
comment
Разве вы не можете заменить последнюю строку на U <- apply(U, 1, cumprod)? - person Joshua Ulrich; 21.03.2013
comment
вы, вероятно, можете добиться большего, используя cumsum на правильно масштабированной rnorm() матрице, а затем возведя в степень - person Ben Bolker; 21.03.2013
comment
@ Ferdinand.kraft какие-либо указания, если и насколько это решение быстрее, чем то, что предлагает OP? - person Paul Hiemstra; 21.03.2013
comment
@PaulHiemstra, примерно на 2 порядка быстрее. См. Правку выше. - person Ferdinand.kraft; 21.03.2013
comment
Всем привет! Большое спасибо за внимание. Использование этого кода действительно намного быстрее. - person PereMkB; 22.03.2013
comment
Сейчас я просто думаю, как я могу использовать эти советы, если цена акций следует более сложному стохастическому процессу. - person PereMkB; 22.03.2013