Variance-Gamma Model

gistfile1.txt

"""
I'm starting to like the Variance-Gamma model quite a bit.
Here is a quick and dirty Octave (Matlab) code for generating the IV skew given VG parameters
I took the parameters from the Madan, Carr and Chang paper.
Instead of evaluating a double integral I use the mixing formula, which results in a single integral.
From Linkedin profile of Frido Rolloos
""


import numpy as np
from scipy.integrate import quad
from scipy.stats import gamma

def vgcall():
    # VG parameters from Madan, Carr and Chang
    sigma = 0.1214
    nu = 0.1686
    theta = -0.1436  # For a perfectly symmetric smile use theta = -0.5*sigma**2

    S0 = 1
    T = 0.5
    r = 0

    ret = []

    for K in np.arange(0.75, 1.26, 0.01):
        def vgcallint(x):
            a = T / nu  # shape parameter
            b = nu  # scale parameter
            omega = np.log(1 - theta * nu - 0.5 * sigma ** 2 * nu) / nu
            xi = np.exp(omega * T + theta * x + 0.5 * sigma ** 2 * x)
            return black_scholes_price(S0 * xi, K, r, T, sigma * np.sqrt(x / T)) * gamma.pdf(x, a, scale=b)

        tmp, _ = quad(vgcallint, 0, np.inf)
        ret.append(black_scholes_implied_volatility(S0, K, r, T, tmp))

    K = np.arange(0.75, 1.26, 0.01)
    import matplotlib.pyplot as plt
    plt.plot(K, ret)
    plt.xlabel("K")
    plt.ylabel("IV")
    plt.show()

def black_scholes_price(S, K, r, T, sigma):
    from scipy.stats import norm
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)

def black_scholes_implied_volatility(S, K, r, T, C):
    from scipy.optimize import fsolve
    
    def f(sigma):
        return C - black_scholes_price(S, K, r, T, sigma)
    
    return fsolve(f, 0.2)[0]

vgcall()