1 
>   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

こんにちは、20000 行を超える上記のデータフレームを R で作成しました。ボラティリティが GARCH(1,1) プロセスに従うと仮定して、オプションの理論価格を計算するコードを書きました。コードは正常に動作しますが、非常に遅いです。それを高速化またはベクトル化する機会があるのだろうか？計算してみましたが、初心者のRユーザーとして失敗しました。計算はモンテカルロシミュレーションで行われます。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

4

2 に答える 2

0

それは私の質問に対するベクトル化された解決策です！

#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

}
)
于 2013-02-25T17:56:56.193 に答える
0

一般に、データ フレームのインデックス作成 (多くの場合) には、多くの時間がかかります。これをベクトル化することを検討してください。特に、これらすべての 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
}
于 2013-02-12T16:12:04.987 に答える