使用蒙特卡洛模拟遵循 Garch 过程的波动性期权定价

Question

>   Name    Date    Close   CP  ttmDaysW    ttm Strike  Fut Wibor   lambda  omega   alpha   beta    sigma
1   OW20C1330   2011-01-19  0.60    c   42  0.1673307   3300    2768    0.0425  0.03985676  1.205098e-06    0.05403404  0.9426635   0.010935144
2   OW20C1330   2011-02-16  0.21    c   22  0.0876494   3300    2703    0.0435  0.03285167  5.852091e-07    0.05208226  0.9462142   0.008209948
3   OW20C2150   2011-12-08  745.65  c   71  0.2828685   1500    2233    0.0499  0.05490974  1.213260e-06    0.06837361  0.9296792   0.018583414
4   OW20C2150   2011-12-09  720.80  c   70  0.2788845   1500    2262    0.0499  0.05119041  1.212956e-06    0.06813476  0.9299286   0.019143222

嗨，我在 R 中创建了上面的数据框，它有超过 20000 行。我编写了一个代码来计算期权的理论价格，假设波动率遵循 GARCH(1,1) 过程。该代码工作正常，但非常缓慢。我想知道天气有没有机会加快速度或矢量化？我试图解决这个问题，但作为一个 R 初学者，我失败了。计算是由 Monte Caro Simulation 完成的。 OW 是我的 Data.Frame

#Monte Carlo Garch(1,1)
nsim=10000
for (i in 1:nrow(OW)){
  iopt<-ifelse(OW$CP[i]=="c",1,-1)
  sum=0
  for (j in 1:nsim){
    Sigma2t<-(OW$sigma[i])^2
    Eps<-rnorm(1)*OW$sigma[i]

    sumSigma2t=0
    sumEps=0

    for (k in 1:OW$ttmDaysW[i]){
      Sigma2t= OW$omega[i] +OW$alpha[i]*(Eps-OW$lambda[i]*sqrt(Sigma2t))^2+OW$beta[i]*Sigma2t
      Eps <- rnorm(1)*sqrt(Sigma2t)

      sumEps=sumEps+Eps
      sumSigma2t = sumSigma2t + Sigma2t

    }
    Ft<-OW$Fut[i]*exp(-0.5*sumSigma2t+sumEps)
    payoff <- max(c(iopt * (Ft - OW$Strike[i]), 0))
    sum<-sum+payoff  
  } 
  OW$G[i] = exp(-OW$Wibor[i] * OW$ttm[i]) * sum / nsim
}

我在我的问题上只找到了这个帮助：Simulation of GARCH in R

Answer 1

这是我问题的矢量化解决方案！

#Monte Carlo Garch(1,1)
N=10000
system.time(for (i in 1:nrow(OW)){
  iopt<-ifelse(OW$CP[i]=="c",1,-1)
  h0<- (OW$sigma[i])^2
  omega <- OW$omega[i]
  alpha1 <- OW$alpha[i]
  lambda <- OW$lambda[i]
  beta1 <- OW$beta[i]
  Fu <- OW$Fut[i]
  Str <- OW$Strike[i]
  horizon <- OW$ttmDaysW[i]
  Wibor<-OW$Wibor[i]
  ttm<-OW$ttm[i]
  ret <- et <- ht <- matrix(NA, nc = horizon, nr = N)
  zt  <- matrix(rnorm(N * horizon, 0, 1), nc = horizon)
  hB<-matrix(rep(h0,N),nr=N)
  eB<-matrix(rnorm(N, 0, 1), nc=1) * sqrt(hB)
  Fut<-matrix(rep(Fu,N),nr=N)
  Strike<-matrix(rep(Str,N),nr=N)
  for(j in 1:horizon){
    ifelse(j>1,
           ht[, j ] <- omega + alpha1 * (et[, j-1]-lambda*sqrt(ht[, j-1])) ^ 2 + beta1 * ht[, j-1],
           ht[, j ] <- omega + alpha1 * (eB -lambda*sqrt(hB))^ 2 + beta1 * hB
    )
    et[, j]  <- zt[, j] * sqrt(ht[, j])
  }
  Ft<-Fut*exp(-0.5*rowSums(ht)+rowSums(et))
  payoff<-pmax(iopt * (Ft - Strike),0)
  OW$G[i] = exp(-Wibor *ttm) * sum(payoff) / N

}
) 

Answer 2

通常，索引数据帧（您经常这样做）会花费大量时间。考虑将其向量化。特别是因为您正在计算所有这些 OW$something[i] 超过 2 亿次（因为所有嵌套的 for 循环），实际上它们只需要被调用 10,000 次（nsim 次）。看看这是否运行得更快。

nsim=10000
for (i in 1:nrow(OW)){
  iopt<-ifelse(OW$CP[i]=="c",1,-1)
  sum=0
  OWSigma <- OW$sigma[i]
  OWOmega <- OW$omega[i]
  OWAlpha <- OW$alpha[i]
  OWLambda <- OW$lambda[i]
  OWBeta <- OW$beta[i]
  OWFut <- OW$Fut[i]
  OWStrike <- OW$Strike[i]
  OWTtmDaysW <- OW$ttmDaysW[i]
  for (j in 1:nsim){
    Sigma2t<-(OWSigma)^2
    Eps<-rnorm(1)*OWSigma

    sumSigma2t=0
    sumEps=0

    for (k in 1:OWTtmDaysW){
      Sigma2t= OWOmega +OWAlpha*(Eps-OWLambda*sqrt(Sigma2t))^2+OWBeta*Sigma2t
      Eps <- rnorm(1)*sqrt(Sigma2t)

      sumEps=sumEps+Eps
      sumSigma2t = sumSigma2t + Sigma2t

    }
    Ft<-OWFutexp(-0.5*sumSigma2t+sumEps)
    payoff <- max(c(iopt * (Ft - OWStrike, 0))
    sum<-sum+payoff  
  } 
  OW$G[i] = exp(-OW$Wibor[i] * OW$ttm[i]) * sum / nsim
}