Calculating implied vol using vanna volga

Question

I'm trying to write an R script which takes in the 25d calls and puts along with ATM and creates a vol. smile.

I've attempted to use the method put forward by Castagna and Mercurio to calculate the vanna-volga volatility. The expression put forward, for a European Call is:

C(K) = CBS(K)+ w(K1)[CMkt(K1)-CBS(K1)]+w(K3)[CMkt(K3)-CBS(K3)]

where w1 and w3 are calculated as per the code below.

BSoption <- function(type, S, X, t, r, rf, v)
{

  d1 <- (log(S/X) + (r-rf + 0.5 * v^2) * t)/(v * sqrt(t))
  d2 <- d1 - v*sqrt(t)
  if (type == "c"){
    pnorm(d1)*S*exp(-rf*t) - pnorm(d2)*X*exp(-r*t)
  } else {
    pnorm(-d2)*X*exp(-r*t) - pnorm(-d1)*S*exp(-rf*t)
  }
}

vega <- function(S, X, t, r, rf, v)
{
  d1 <- (log(S/X) + (r -rf + 0.5 * v^2) * t)/(v * sqrt(t))
  Np <- (exp(-d1^2/2))/ sqrt(2 * pi) 
  (S * exp(-rf*t)*sqrt(t) * Np)/100
}

implied.vol <-
  function(type, S, X, t, r, rf, market){
    sig <- 0.20
    sig.up <- 1
    sig.down <- 0.001
    count <- 0
    err <- BSoption(type, S, X, t, r, rf, sig) - market 

    ## repeat until error is sufficiently small or counter hits 1000
    while(abs(err) > 0.00001 && count<1000){
      if(err < 0){
        sig.down <- sig
        sig <- (sig.up + sig)/2
      }else{
        sig.up <- sig
        sig <- (sig.down + sig)/2
      }
      err <- BSoption(type, S, X, t, r, rf, sig) - market
      count <- count + 1
    }

    ## return NA if counter hit 1000
    if(count==1000){
      return(NA)
    }else{
      return(sig)
    }
  }

S <- 0.906
X <- seq(0.7,1.2,0.01)
t <- 1
r <- 0.0507
rf <- 0.047#-log(0.9945049)/t
ATMcost <- BSoption("c",S,XATM,t,r,rf,vATM)

v25p <- 0.13575#vv.inputs$Vol[vv.inputs$Skew == -0.25]
vATM <- 0.132#vv.inputs$Vol[vv.inputs$Skew == 0.0]
v25c <- 0.13425#vv.inputs$Vol[vv.inputs$Skew == 0.25]

X25p <- 0.8350575#BSStrikeFromDelta("p",S,v25p,t,r,rf,0.25)
XATM <- S
X25c <- 1.000846#BSStrikeFromDelta("c",S,t,r,v25c,rf,0.25)

w1 <- (vega(S,X,t,r,rf,vATM)/vega(S,X25p,t,r,rf,v25p))*((log(XATM/X)*log(X25c/X))/(log(XATM/X25p)*log(X25c/X25p)))
w2 <- (vega(S,X,t,r,rf,vATM)/vega(S,XATM,t,r,rf,vATM))*((log(X/X25p)*log(X25c/X))/(log(XATM/X25p)*log(X25c/XATM)))
w3 <- (vega(S,X,t,r,rf,vATM)/vega(S,X25c,t,r,rf,v25c))*((log(X/X25p)*log(X/XATM))/(log(X25c/X25p)*log(X25c/XATM)))

VV.price <- ATMcost + w1*(BSoption("c",S,X25p,t,r,rf,v25p)-BSoption("c",S,X25p,t,r,rf,vATM)) + w3*(BSoption("c",S,X25c,t,r,rf,v25c)-BSoption("c",S,X25c,t,r,rf,vATM))

VV.vol <- 0
for(i in 1:length(X)){
  VV.vol[i] <- implied.vol("c",S,X[i],t,r,rf,VV.price[i])
}
plot(X,VV.price)
plot(X,VV.vol)

As you can see, I calculate the vanna-volga price of the option, which is consistent with the given quotes but, when I plot the volatility 'smile', it actually is not a smile at all and I'm not sure where I'm going wrong.

Any help is much appreciated.

Answer 1

In the line below all I had to do was replace X[i] with XATM because the vanna volga method creates a price for an option that satisfies the hedged vega,vanna,volga ATM so you're always getting vol. with ATM strike.

VV.vol[i] <- implied.vol("c",S,X[i],t,r,rf,VV.price[i])

Answer 2

Actually the best way to answer this question was to use the original Mercurio and Castagna paper which has a proof and explanation for the best way to calculate implied vol.

See paper here:

http://www.fabiomercurio.it/consistentfxsmile.pdf