Model-Free Dirac Delta Integration for Implied Volatility.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model-Free Dirac Delta Integration for Implied Volatility\n",
    "\n",
    "Author: Zhenyu Cui, Justin Kirkbyb, Duy Nguyenc and Stephen Taylord\n",
    "\n",
    "DOI: [10.13140/RG.2.2.21031.29602](http://dx.doi.org/10.13140/RG.2.2.21031.29602)\n",
    "\n",
    "From the paper:\n",
    "\n",
    "\"The method is based on the novel use of the Dirac Delta function, corresponding delta sequences, and the change of variable technique. The formula is expressed through either a limit or as an infinite series of elementary functions, and we establish that the proposed formula converges to the true implied volatility value. In numerical experiments , we verify the convergence of the formula, and consider several benchmark cases, for which the data generating processes are respectively the stochastic volatility inspired (SVI) model, the stochastic alpha beta rho (SABR) model. We also establish an explicit formula for the implied volatility expressed directly in terms of respective model parameters , and use the Heston model to illustrate this idea. The delta sequence and change of variable technique that we develop is of independent interest and can be used to solve inverse problems arising in other applications.\"\n",
    "\n",
    "## Formula\n",
    "\\begin{equation}\n",
    "    \\sigma_{BS}(C) = \\lim_{\\epsilon \\rightarrow 0} \\frac{S_0 \\sqrt{T}}{2 \\pi \\sqrt{2 \\epsilon}} \\cdot \\int_{\\sigma_l}^{\\sigma_u} s \\cdot \\exp \\left( -\\frac{\\left(f_{BS}(s)-C\\right)^2}{4\\epsilon} -\\frac{d_1^2(s)}{2}\\right) ds\\\\\n",
    "\\end{equation}\n",
    "where $\\sigma_l$ and $\\sigma_u$ are the upper bound and lower bound that contains the true $\\sigma_{BS}$ which can be narrow down using the following formulas (Tehranchi 2016.Uniform bounds for Black-Scholes implied volatility. SIAM Journal on Financial Mathematics)\n",
    "\\begin{equation}\n",
    "\\sigma_u = \\frac{-2}{\\sqrt{T}}\\mathcal{N}^{-1}\\left(\\frac{1-c}{1+e^k}\\right)\\\\\n",
    "\\sigma_l = \\begin{cases}\n",
    "\\frac{-2}{\\sqrt{T}}\\mathcal{N}^{-1}\\left(\\frac{1-c}{2e^k}\\right)\\quad if\\space k<0\\\\\n",
    "\\frac{-2}{\\sqrt{T}}\\mathcal{N}^{-1}\\left(\\frac{1-c}{2}\\right)\\quad if\\space k\\geq 0\n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "Note that this formula is for call option only, but can be used for put options by put-call parity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "import scipy.integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BS formulas\n",
    "def d1_fun(S,K,T,r,q,sigma):\n",
    "    return (np.log(S/K)+(r+0.5*sigma**2)*T)/sigma/np.sqrt(T)\n",
    "\n",
    "def d2_fun(S,K,T,r,q,sigma):\n",
    "    return (np.log(S/K)+(r-0.5*sigma**2)*T)/sigma/np.sqrt(T)\n",
    "\n",
    "def black_scholes_call(S,K,T,r,q,sigma):\n",
    "    forward = S*np.exp((r-q)*T)\n",
    "    discount_factor = np.exp(-r*T)\n",
    "    d1 = d1_fun(S,K,T,r,q,sigma)\n",
    "    d2 = d1-sigma*np.sqrt(T)\n",
    "    N = lambda x: sp.stats.norm.cdf(x)\n",
    "    return np.exp(-q*T)*S*N(d1)-np.exp(-r*T)*K*N(d2)\n",
    "\n",
    "def black_scholes_put(S,K,T,r,q,sigma):\n",
    "    forward = S*np.exp((r-q)*T)\n",
    "    discount_factor = np.exp(-r*T)\n",
    "    d1 = d1_fun(S,K,T,r,q,sigma)\n",
    "    d2 = d1-sigma*np.sqrt(T)\n",
    "    N = lambda x: sp.stats.norm.cdf(x)\n",
    "    return -np.exp(-q*T)*S*N(-d1)+np.exp(-r*T)*K*N(-d2)\n",
    "\n",
    "# IV Boundary\n",
    "def c_fun(C,S,T,q):\n",
    "    return C/(S*np.exp(-q*T))\n",
    "\n",
    "def k_fun(S,K,T,r,q):\n",
    "    return np.log(K/(S*np.exp((r-q)*T)))\n",
    "\n",
    "def sigma_bounds_call(C,S,K,T,r,q):\n",
    "    \"\"\"The bounds for call options\"\"\"\n",
    "    c = c_fun(C,S,T,q)\n",
    "    k = k_fun(S,K,T,r,q)\n",
    "    invN = lambda x: sp.stats.norm.ppf(x)\n",
    "    sigma_upper = -2/np.sqrt(T)*invN((1-c)/(1+np.exp(k)))\n",
    "    if k<0:\n",
    "        sigma_lower = -2/np.sqrt(T)*invN((1-c)/(2*np.exp(k)))\n",
    "    else:\n",
    "        sigma_lower = -2/np.sqrt(T)*invN((1-c)/(2))\n",
    "    return sigma_lower,sigma_upper\n",
    "\n",
    "def sigma_bounds_put(P,S,K,T,r,q):\n",
    "    \"\"\"The bounds for put options\"\"\"\n",
    "    C = P + S - np.exp(-r*T)*K\n",
    "    c = c_fun(C,S,T,q)\n",
    "    k = k_fun(S,K,T,r,q)\n",
    "    invN = lambda x: sp.stats.norm.ppf(x)\n",
    "    sigma_upper = -2/np.sqrt(T)*invN((1-c)/(1+np.exp(k)))\n",
    "    if k<0:\n",
    "        print((1-c)/(2*np.exp(k)))\n",
    "        sigma_lower = -2/np.sqrt(T)*invN((1-c)/(2*np.exp(k)))\n",
    "    else:\n",
    "        sigma_lower = -2/np.sqrt(T)*invN((1-c)/(2))\n",
    "    return sigma_lower,sigma_upper\n",
    "\n",
    "def dirac_delta_implied_vol(C,S,K,T,r,q,sigma_lower,sigma_upper,epsilon):\n",
    "    def integral_fun(s):\n",
    "        return s*np.exp(-(black_scholes_call(S,K,T,r,q,s)-C)**2/4/epsilon-(d1_fun(S,K,T,r,q,s)**2/2))\n",
    "    coef = S*np.sqrt(T)/(2*np.pi*np.sqrt(2*epsilon))\n",
    "    integral,err = sp.integrate.quad(integral_fun,sigma_lower,sigma_upper)\n",
    "    return coef*integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "K_range = np.array([100 + 10*k for k in range(11)])\n",
    "T_range = np.array([1/48,1/12,3/12,6/12,1])\n",
    "r_range = np.array([-0.01,0,0.01,0.05])\n",
    "q_range = np.array([0,0.01,0.05,0.1])\n",
    "sigma_range = np.array([0.1,0.3,0.5,0.7])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Have a Try\n",
    "Run the following cell, it would test a random market case each time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "###################################\n",
      "Static Parameters\n",
      "###################################\n",
      "sigma: 0.5\n",
      "r:     -0.01\n",
      "q:     0.01\n",
      "K:     200\n",
      "T:     0.5\n",
      "###################################\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=0.4\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "S:     80.00\n",
      "call:  0.06\n",
      "put:   121.46\n",
      "lower_vol: 0.0026972877882851733\n",
      "upper_vol: 1.6195221025849988\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.3473\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.3688\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.4732\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.5193\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=0.7\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "S:     140.00\n",
      "call:  4.54\n",
      "put:   66.24\n",
      "lower_vol: 0.11564630045059643\n",
      "upper_vol: 0.7459388464124055\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.5083\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.5086\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.5086\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.0000\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=0.8\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "S:     160.00\n",
      "call:  9.59\n",
      "put:   51.39\n",
      "lower_vol: 0.2137042835372123\n",
      "upper_vol: 0.6046741171437177\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.5058\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.5059\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.0000\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.0000\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=0.9\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "S:     180.00\n",
      "call:  17.11\n",
      "put:   39.01\n",
      "lower_vol: 0.3393732302145045\n",
      "upper_vol: 0.5263205077351348\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.5026\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.5035\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.5035\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.5035\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=1\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "S:     200.00\n",
      "call:  27.07\n",
      "put:   29.07\n",
      "lower_vol: 0.48460298264555607\n",
      "upper_vol: 0.5001499950522819\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.2377\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.2614\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.2874\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.3619\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=1.1\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "0.44415977768728554\n",
      "S:     220.00\n",
      "call:  39.25\n",
      "put:   21.35\n",
      "lower_vol: 0.3776647618050493\n",
      "upper_vol: 0.5144326346408324\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.4819\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.4994\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.4994\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.4994\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=1.2\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "0.4584683596231978\n",
      "S:     240.00\n",
      "call:  53.30\n",
      "put:   15.50\n",
      "lower_vol: 0.27376799067298074\n",
      "upper_vol: 0.5581274057115231\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.4977\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.4977\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.4977\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.0000\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=1.3\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "0.469036550663551\n",
      "S:     260.00\n",
      "call:  68.85\n",
      "put:   11.15\n",
      "lower_vol: 0.19681409680450482\n",
      "upper_vol: 0.6218915282991315\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.4961\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.4961\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.4961\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.4961\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "Moneyness S/K=1.6\n",
      "$$$$$$$$$$$$$$$$$$$$$\n",
      "0.4860359428203269\n",
      "S:     320.00\n",
      "call:  121.42\n",
      "put:   4.02\n",
      "lower_vol: 0.0708626081809273\n",
      "upper_vol: 0.8696035199720772\n",
      "epsilon: 0.1      , vol: 0.50, dirac_delta_vol 0.4917\n",
      "epsilon: 0.01     , vol: 0.50, dirac_delta_vol 0.4921\n",
      "epsilon: 0.001    , vol: 0.50, dirac_delta_vol 0.4921\n",
      "epsilon: 0.0001   , vol: 0.50, dirac_delta_vol 0.0000\n"
     ]
    }
   ],
   "source": [
    "K = np.random.choice(K_range)\n",
    "T = np.random.choice(T_range)\n",
    "r = np.random.choice(r_range)\n",
    "q = np.random.choice(q_range)\n",
    "sigma = np.random.choice(sigma_range)\n",
    "\n",
    "print('###################################')\n",
    "print('Static Parameters')\n",
    "print('###################################')\n",
    "print('sigma: {}'.format(sigma))\n",
    "print('r:     {}'.format(r))\n",
    "print('q:     {}'.format(q))\n",
    "print('K:     {}'.format(K))\n",
    "print('T:     {}'.format(T))\n",
    "print('###################################')\n",
    "for moneyness in [0.4,0.7,0.8,0.9,1,1.1,1.2,1.3,1.6]:\n",
    "    print('$$$$$$$$$$$$$$$$$$$$$')\n",
    "    print('Moneyness S/K={}'.format(moneyness))\n",
    "    print('$$$$$$$$$$$$$$$$$$$$$')\n",
    "    S = K*moneyness\n",
    "    call_price = black_scholes_call(S,K,T,r,q,sigma)\n",
    "    put_price = black_scholes_put(S,K,T,r,q,sigma)\n",
    "    call_sigma_bounds = sigma_bounds_call(call_price,S,K,T,r,q)\n",
    "    put_sigma_bounds = sigma_bounds_put(put_price,S,K,T,r,q)\n",
    "    print('S:     {:.2f}'.format(S))\n",
    "    print('call:  {:.2f}'.format(call_price))\n",
    "    print('put:   {:.2f}'.format(put_price))\n",
    "    print('lower_vol: {}'.format(call_sigma_bounds[0]))\n",
    "    print('upper_vol: {}'.format(call_sigma_bounds[1]))\n",
    "    for epsilon in [0.1,0.01,0.001,0.0001]:\n",
    "        dirac_delta_vol = dirac_delta_implied_vol(call_price,S,K,T,r,q,call_sigma_bounds[0],call_sigma_bounds[1],epsilon)\n",
    "        print('epsilon: {:<9}, vol: {:.2f}, dirac_delta_vol {:.4f}'.format(epsilon,sigma,dirac_delta_vol))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "Simple in form and good to know a close-form model-free formula on implied volatility. Indeed, getting the implied volatility from an arbitrary parameterized model is valuable especially in deriative's trading. Investment banks tend to use their own models in valuation but ultimately would be interpreted as implied vol and compare it market bid/offer from rivals. So this framework should be a long-last and forever topic in Financial Engineering.\n",
    "\n",
    "### On clarity:\n",
    "The implied volatility barrier formula (Lemma 1) is only for **call option case**. The put case boundaries can be addressed easily by call-put parity but better to claim it explicitly\n",
    "\n",
    "### On boundary estimation\n",
    "The formula is simple (easy to implement) but not very tight, do we have better estimates on that?\n",
    "\n",
    "### On numerical implementation\n",
    "\n",
    "Because of epsilon and Dirac Delta function, it requires high precision in integration (finer steps), otherwise would not be easy to get stable solutions. See following cases\n",
    "1. early expiry\n",
    "2. too small epsilon would return zero without using sophisticated numerical integration (step should adaptive to epsilon)\n",
    "3. deep ITM or OTM\n",
    "Therefore, may be practical to elaborate how to numerically solve the integration part."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}