Option Pricing with volatility following a Garch process by use of Monte-Carlo Simulations

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

Hi, I created the above dataframe in R which has above 20000 rows. I wrote a code to compute theoretical prices of Options assuming that volatility follow a GARCH(1,1) process. The code works fine but is VERY sluggish. I wonder weather there is any chance to speed it up or Vectorize? I've tried to work it out, but I failed as a beginning R user.Computation is done by Monte Caro Simulation. OW is my 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
}

I have found only this help on my question:Simulation of GARCH in R

score 0 · Answer 1 · answered Feb 12 '13 at 16:12

Generally, indexing data frames (which you do a lot), takes a lot of time. Consider vectorizing this. Especially since you are computing all those OW$something[i]'s more than 200 million times (since all the nested for loops) where they actually only have to be called 10,000 times (nsim times). See if this runs any quicker.

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
}

Thanks. It speeded up the code quite well. As you said, i'm doing my best at vectorizing this. — krson, Feb 13 '13 at 10:06
then you can also click "accept this answer" if this solved your problem ;-) — Gx1sptDTDa, Feb 13 '13 at 13:50
i would accept but i stil expect a vectorized solution as the one that is satisfactory to me, because this takes over a week to get my results ( i have couple similar models to compute). Bus still very thankful for what you have done to me :) — krson, Feb 15 '13 at 12:32

score 0 · Accepted Answer · answered Feb 25 '13 at 17:56

It's vectorized solution to my question!

#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

}
)

Option Pricing with volatility following a Garch process by use of Monte-Carlo Simulations

2 Answers2