Моделирование броуновского движения с использованием R - PullRequest
0 голосов
/ 09 мая 2020

Моделирование броуновского движения в инвертированном времени [0,100], и пути были нарисованы путем моделирования n = 1000 точек. Я генерирую следующий код:

 n <- 1000
 t <- 100
 bm <- c(0, cumsum(rnorm(n,0,sqrt(t/n))))
 steps <- seq(0,t,length=n+1)
 plot(steps,bm,type="l")

Как я могу смоделировать 50 примеров траекторий стандартного броуновского движения и показать каждый путь другим цветом, как набор траекторий?

Я думаю это будет что-то вроде replicate(50,bm), но когда я это сделаю, в xy.coords будет ошибка. Спасибо за помощь!

Моделирование Броуновского моста на [0,1], и пути были нарисованы путем моделирования n = 1000 точек. Я генерирую следующий код

n <- 1000
t <- seq(0,1,length=n)
No.Ex<-10
bm <- c(0,cumsum(rnorm(n-1,0,1)))/sqrt(n)
B = replicate(No.Ex,{
  bb <- bm - t*bm[n]
})
matplot(B, type = "l", col = cols, lty = 1)

Код для генерации примеров путей для Geometri c Brownian Motion

simGBM<- function(P0, mu, sigma, T, nSteps, nRepl){
  dt<- T/nSteps
  muT<- (mu-sigma^2/2)*dt
  sigmaT<- sqrt(dt)*sigma
  pathMatrix<- matrix(nrow = nRepl, ncol = nSteps+1)
  pathMatrix[,1]<- P0
  for(i in 1:nRepl){
    for(j in 2:(nSteps+1)){
      pathMatrix[i,j]<- pathMatrix[i,j-1]*exp(rnorm(1, muT, sigmaT))
    }
  }
  return(pathMatrix)
}

P0<- 1 #initial price
mu<- 0.1 #drift
sigma<- 0.5 #volatility
T<- 100/360 #100 days of a commercial year
nSteps<- 50 #No of steps
nRepl<- 100 #No of replications

paths<- simGBM(P0, mu, sigma, T, nSteps, nRepl)
yBounds<- c(min(paths),max(paths)) #bounds of simulated prices

plot(paths[1,], ylim = yBounds, type = 'l',col = 1, main = "Simulation of sample paths of GBM", xlab = "Time", ylab = "Price")
for(k in 2:numRepl) lines(paths[k,], col = k)

Я пытаюсь использовать функцию matplot, но не могу создать то же самое график

cols = rainbow(nSteps)
matplot(paths, ylim = yBounds, type = "l", col = cols, lty = 1, main = "Simulation of sample paths of GBM", xlab = "Time", ylab = "Price")

1 Ответ

2 голосов
/ 09 мая 2020

Как насчет этого

n = 1000
t = 100
No.Ex = 10
steps = seq(0,t,length=n+1)
A = replicate(No.Ex, {
  bm <- c(0, cumsum(rnorm(n,0,sqrt(t/n))))
}) 

cols = rainbow(No.Ex)
matplot(A, type = "l", col = cols, lty = 1)

Я изменил свой ответ и включил предложение Стефана Лорана matplot. Это дает следующее изображение.

enter image description here

EDIT:

Чтобы ответить на ваш вопрос в комментариях, я думаю, вам следует сохранить мой исходный код для bm что составляет bm <- c(0, cumsum(rnorm(n,0,sqrt(t/n)))). Тогда все работает очень хорошо! Спасибо, что указали на красивую matplot команду @Stephane Laurent.

EDIT2: Я только что понял, что вы задали новый вопрос в отношении Коричневого моста. Вы можете попробовать этот код

n <- 1000
t <- seq(0,1,length=n)
No.Ex<-10
B = replicate(No.Ex,{
  bm <- c(0, cumsum(rnorm(n - 1,0,sqrt(t/n))))
  bb <- bm - t*rep(bm[length(bm)], length.out = length(bm))
})
matplot(B, type = "l", col = cols, lty = 1)

Это дает

enter image description here

Кроме того, для Geometri c Brownian Motian попробуйте эту модификацию ваш код с меньшим количеством репликаций

simGBM<- function(P0, mu, sigma, T, nSteps, nRepl){
  dt<- T/nSteps
  muT<- (mu-sigma^2/2)*dt
  sigmaT<- sqrt(dt)*sigma
  pathMatrix<- matrix(nrow = nRepl, ncol = nSteps+1)
  pathMatrix[,1]<- P0
  for(i in 1:nRepl){
    for(j in 2:(nSteps+1)){
      pathMatrix[i,j]<- pathMatrix[i,j-1]*exp(rnorm(1, muT, sigmaT))
    }
  }
  return(pathMatrix)
}

P0<- 1 #initial price
mu<- 0.1 #drift
sigma<- 0.5 #volatility
T<- 100/360 #100 days of a commercial year
nSteps<- 50 #No of steps
nRepl<- 10 #No of replications

paths<- simGBM(P0, mu, sigma, T, nSteps, nRepl)
yBounds<- c(min(paths),max(paths)) #bounds of simulated prices

plot(paths[1,], ylim = yBounds, type = 'l',col = 1, main = "Simulation of sample paths of GBM", xlab = "Time", ylab = "Price")
for(k in 2:nRepl) lines(paths[k,], col = k)

cols = rainbow(nSteps)
matplot(paths, ylim = yBounds, type = "l", col = cols, lty = 1, main = "Simulation of sample paths of GBM", xlab = "Time", ylab = "Price")

На моей машине это дает

enter image description here

...