# Calculating implied vol using vanna volga

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.

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