在R中模拟一个过程n次答案

【问题标题】：Simulating a process n times in R在R中模拟一个过程n次
【发布时间】：2021-05-22 12:31:56
【问题描述】：

我编写了一个 R 脚本（来自 here）模拟股票价格的几何布朗运动的路径，我需要模拟运行 1000 次，以便生成 1000 条过程路径 Ut = Ste^-mu*t，通过离散化从 Ut 导出的运动定律，这是解决 here 问题的底线。

该过程还具有 n = 252 步骤和离散化步骤 = 1/252，还有 sigma = 0.4 的风险和瞬时漂移 mu，我将其视为零，尽管我对此不确定。我正在努力模拟流程的 1000 条路径，但能够生成一条路径，我不确定我需要更改哪些变量，或者我的 for 循环中是否存在限制我生成所有 1000 条路径的问题。是否也可能是脚本在模拟 252 实现的每个单独点，而不是模拟整个过程？如果是这样，这会限制我生成所有 1000 条路径吗？我生成的定义为 U 的数组是否也有可能没有被我正确生成？ U[0] 必须等于 1，第一个实现 U(1) = 1 也必须等于 1。代码如下，我很想弄清楚这一点，因此不胜感激。

#Simulating  Geometric Brownian motion (GMB)
tau <- 1 #time to expiry
N <- 253 #number of sub intervals
dt <- tau/N #length of each time sub interval
time <- seq(from=0, to=N, by=dt) #time moments in which we simulate the process
length(time) #it should be N+1

mu <- 0 #GBM parameter 1
sigma <- 0.4 #GBM parameter 2
s0 <- 1 #GBM parameter 3

#simulate Geometric Brownian motion path
dwt <- rnorm(N, mean = 0, sd = 1) #standard normal sample of N elements
dW <- dwt*sqrt(dt) #Brownian motion increments
W <- c(0, cumsum(dW)) #Brownian motion at each time instant N+1 elements


#Define U Array and set initial values of U
U <- array(0, c(N,1)) #array of U
U[0] = 1
U[1] <- s0 #first element of U is s0. with the for loop we find the other N elements

for(i in 2:length(U)){
  U[i] <- (U[1]*exp(mu - 0.5*sigma^2*i*dt + sigma*W[i-1]))*exp(-mu*i)
}

#Plot 
plot(ts(U), main = expression(paste("Simulation of Ut")))

【问题讨论】：

标签： r simulation stochastic-process

【解决方案1】：

这个问题很难回答，因为有很多不清楚的地方，至少对我来说是这样。

首先，length(time) 等于 64010，而不是 N + 1，后者将是 254。

如果我理解正确，布朗运动函数会返回给定时间的一维位置。因此，每次计算这个位置就足够了：

s0*exp((mu - 0.5*sigma^2)*time + sigma*rnorm(length(time),0,time))

但是，这计算的是 64010 点，而不是 253。如果你复制它 1000 次，它会给出 64010000 点，这是相当多的。

> B <- 1000
> res <- replicate(B, {
+   s0*exp((mu - 0.5*sigma^2)*time + sigma*rnorm(length(time),0,time))
+ })
> length(res)
[1] 64010000
> dim(res)
[1] 64010  1000

我知道我错过了第二部分，即解释 here 的部分，但我实际上并不完全理解您在那里需要什么。如果你能画出公式，也许我可以帮助你。

一般来说，避免在 R 中使用for 循环来迭代向量。 R 是一种矢量化语言，没有必要这样做。如果你想运行相同的代码B 次，replicate(B,{ your code }) 函数就是你的朋友。

【讨论】：