Volatility Surface Modeling for Equity markets


i am trying to construct volatility surface projection for stocks via the Vanna-Volga Implied Volatility method.

The Vol Smile i am getting is way wrong if i go by 3 inputs for the given vol curve:

Input TYPESigmaStrike
25 Delta PUT0.246814650
25 Delta Call0.208814950

Screenshot 2021-03-23 at 4.29.40 PM.png

Actual Vol Smile: 1.95 Days to expiry, Atm Price is 14850
Screenshot 2021-03-23 at 6.21.13 PM.png

Generated Vol Smile via Vanna-Volga Implied Volatility

I believe since the Volsmile is shifted to the right, this is causing problems, what adjustments should I make to get it done right

Here is the code i am using:

import numpy as np
import scipy.stats as sct
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def volsurface(T,DP,DC,ATM,st,cmp,delta):
    # T is is time in days
    # DP,DC,ATM are Sigma 25 Delta Put,25 Delta Call, And Sigma of ATM
    # st is strike array
    # CMP is ATM
     #delta is =0.25
    T = np.array(T)
    sigma25DeltaPut = np.array(DP)
    sigmaATM = np.array(ATM)
    sigma25DeltaCall = np.array(DC)

    #discount factor Applicable in currencies
    r = np.array([0.00])
    q = np.array([0.00])
    mu = r - q
    Dd = np.exp(-np.multiply(r,T))
    Df = np.exp(-np.multiply(q,T))

    S0 = cmp

    F = S0*np.exp(np.multiply(mu,T))
    alpha = -sct.norm.ppf(delta*np.reciprocal(Df))
    K25DPut = F*np.exp(-(alpha*sigma25DeltaPut*np.sqrt(T)) +0.5*(sigma25DeltaPut**2)*T)
    KATM = F*np.exp(0.5*(sigmaATM**2)*T)
    K25DCall = F*np.exp((alpha*sigma25DeltaCall*np.sqrt(T))+0.5*(sigma25DeltaCall**2)*T)
    def d1(F,K,sigma,t):
        dplus = (np.log(F/K)+0.5*(sigma**2)*t)/(sigma*np.sqrt(t))
        return dplus

    def d2(F,K,sigma,t):
        dminus = d1(F,K,sigma,t) - sigma*np.sqrt(t)
        return dminus

    def VannavolgaImpliedVol(F,K,t,K1,K2,K3,sigma1,sigma2,sigma3):
        #sigma1,2,3 are 25d put,atm 25d call
        y1 = np.array(np.log(K2/K)*np.log(K3/K))/np.array(np.log(K2/K1)*np.log(K3/K1))
        y2 = np.array(np.log(K/K1)*np.log(K3/K))/np.array(np.log(K2/K1)*np.log(K3/K2))
        y3 = np.array(np.log(K/K1)*np.log(K/K2))/np.array(np.log(K3/K1)*np.log(K3/K2))
        P = (y1*sigma1 + y2*sigma2 + y3*sigma3) - sigma2
        Q = y1 * d1(F,K1,sigma1,t) * d2(F,K1,sigma1,t) * ((sigma1-sigma2)**2) \
            + y2 * d1(F,K2,sigma2,t) * d2(F,K2,sigma2,t) * ((sigma2-sigma2)**2) \
            + y3 * d1(F,K3,sigma3,t) * d2(F,K3,sigma3,t) * ((sigma3-sigma2)**2)
        d1d2 = d1(F,K,sigma2,t) * d2(F,K,sigma2,t)
        sigma = sigma2 + (-sigma2 + np.sqrt(sigma2**2 + d1d2 * (2*sigma2*P+Q)))/(d1d2)
        return sigma
    VVImpliedVol = VannavolgaImpliedVol(F, st.reshape(-1,1), T, K25DPut, KATM, K25DCall, sigma25DeltaPut, sigmaATM,sigma25DeltaCall)
    return VVImpliedVol
Last edited:
My 2 cents:
the major problem is following:
your quotes are only using 3 strike points, and they concentrate on the very small strike range [14650, 14950] compared to the big strike range of the actual smile.​
these 3 strike points have a decreasing trend against strikes​
these 3 strike points' vol don't imply that big enough curvature shape like the Actual Smile should have.​
These 3 factors together lead to the declining-slope none-curvature smile shape returned from your model.

To summarize, I don't think using 3 strike points are enough, and using Vanna-Volga Implied Volatility method (this model can only takes 3 points as input) is a not good choice for modelling.
I would suggest you to:
(1) increase the quoted strike points. If you are modeling the listed equity futures market, it should have many more quoted points;
(2) as long as you have many points, even the simplest method like Linear Interpolation should return a similar shape as the Actual Smile, and then you can explore other more complex methods like SABR, SVI, etc.