py/HW3/TestCode_SABR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# SABR Model Test Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import imp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "#import sys\n",
    "#sys.path.insert(sys.path.index('')+1, 'D:/Github/PyFENG')\n",
    "import pyfeng as pf\n",
    "import option_models as opt\n",
    "\n",
    "np.set_printoptions(precision=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Pricing using Hagan's formula (provided)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (1) SABR for $0<\\beta\\le 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "strike = np.linspace(75,125,num=25)\n",
    "forward = 100\n",
    "sigma = 0.2\n",
    "texp = 1\n",
    "vov = 0.5\n",
    "rho = 0.25\n",
    "beta = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create model\n",
    "sabr_bsm = pf.SabrHagan2002(sigma, vov=vov, rho=rho, beta=beta)\n",
    "#sabr_bsm.__dict__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([25.631, 23.751, 21.924, 20.159, 18.465, 16.85 , 15.32 , 13.884,\n",
       "        12.543, 11.302, 10.159,  9.115,  8.165,  7.306,  6.532,  5.838,\n",
       "         5.217,  4.663,  4.17 ,  3.731,  3.342,  2.996,  2.689,  2.417,\n",
       "         2.175]),\n",
       " array([0.204, 0.203, 0.202, 0.201, 0.2  , 0.2  , 0.2  , 0.2  , 0.201,\n",
       "        0.202, 0.203, 0.204, 0.205, 0.206, 0.208, 0.21 , 0.211, 0.213,\n",
       "        0.215, 0.217, 0.219, 0.221, 0.223, 0.225, 0.227]))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is how you price. you can get either price or bsm_vol\n",
    "price = sabr_bsm.price(strike, forward, texp)\n",
    "bsm_vol = sabr_bsm.vol_smile(strike, forward, texp)\n",
    "\n",
    "price, bsm_vol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfHUlEQVR4nO3dfbRddX3n8fcn4amBUOUht2keYRmpGSfgJAWp1iQgGqga0VqJdyIzEiOzJlWpnSGatajUyVqMg2bUguklZGDkDoytWGMN8tRcqBNjCS6aJsSEFJJLIE2UyQgXNJLkO3/sfeFwPefcvU/OPo+f11pnnbMf7/fHDed7f3v/9veniMDMzCyrMc0OwMzM2osTh5mZ5eLEYWZmuThxmJlZLk4cZmaWy3HNDqARzjjjjJg+fXqzw8jtxRdf5OSTT252GA3Tbe0Ft7lbtGubH3300Z9FxJkj13dF4pg+fTqbN29udhi5DQwMMG/evGaH0TDd1l5wm7tFu7ZZ0p5y632pyszMcnHiMDOzXJw4zMwsFycOMzPLxYnDzMxyceIwM+tE/f0wfTqMGZO89/fX7dRdMRzXzKyr9PfD0qXw0kvJ8p49yTJAb+8xn949DjOzTrNixatJY9hLLyXr68CJw8ys0wwO5lufkxOHmVmnmTo13/qcnDjMzDrNypUwbtxr140bl6yvAycOM7NO09sLfX0wbRpIyXtfX11ujINHVZmZdabe3rolipHc4zAzs1ycOMzMLBcnDjMzy8WJw8zMcnHiMDOzXJw4zMwsFycOMzPLpdDEIWmBpB2SdklaXmZ7r6Qt6WujpHPT9VMkbZC0XdI2SZ8qOebzkp6R9Fj6uqzINpiZNV2BJdJrUdgDgJLGAjcBlwB7gUckrYuIx0t2ewqYGxEHJV0K9AEXAIeBz0TEjyWNBx6VdH/Jsasi4saiYjczaxkFl0ivRZE9jvOBXRHxZET8CrgLWFi6Q0RsjIiD6eImYHK6fl9E/Dj9/AKwHZhUYKxmZq2p4BLptSiy5Mgk4OmS5b0kvYlKrgLuGblS0nTgLcCPSlYvk/RRYDNJz+RgmeOWAksBenp6GBgYyBl+8w0NDbVl3LXqtvaC29wtjqXNcwcHUZn1MTjIQ8367xgRhbyADwFrSpYXA1+rsO98kl7F6SPWnwI8CnygZF0PMJakt7QSWDtaLLNnz452tGHDhmaH0FDd1t4It7lbHFObp02LgF9/TZtWp+gqAzZHme/UIi9V7QWmlCxPBp4duZOkWcAaYGFEPFey/njgW0B/RNw9vD4i9kfEkYg4CtxCcknMzKwzFVwivRZFJo5HgBmSzpJ0AnAFsK50B0lTgbuBxRGxs2S9gFuB7RHx5RHHTCxZvBzYWlD8ZmbNV3CJ9FoUdo8jIg5LWgbcS3JpaW1EbJN0dbp9NXAdcDpwc5IrOBwRc4C3kVza+idJj6Wn/FxErAe+KOk8IIDdwCeKaoOZWUsosER6LQqdjyP9ol8/Yt3qks9LgCVljvsBlL0fREQsrnOYZmaWg58cNzOzXJw4zMwsFycOMzPLxYnDzMxyceIwM7NcnDjMzCwXJw4zM8vFicPMrJFabG6NWhT6AKCZmZVowbk1auEeh5lZo7Tg3Bq1cOIwM2uUwcF861uUE4eZWaNMnZpvfYty4jAza5QWnFujFk4cZmaN0oJza9TCo6rMzBqpxebWqIV7HGZmlosTh5mZ5eLEYWZmueROHJIukXR/EcGYmVnrq5g4JF0kaaekIUl3SJopaTNwA/D1xoVoZmatpFqP40vAUuB04K+BTcA3ImJ2RNzdiODMzKz1VEscEREDEXEoIv4G+GlEfKVBcZmZtb4OqHRbi2rPcbxO0gdKllW67F6HmXW1Dql0W4tqieMh4L0VlgNw4jCz7lWt0m23Jo6I+PeNDMTMrK10SKXbWlQdjitprqRZ6ec/kvQXkq6RdGJjwjMza1EdUum2FtWG494E/BfgVkl3AB8BtgJvAdY2JjwzsxbVIZVua1HtHsf8iJgp6STgGWBCRByR9JfAlsaEZ2bWoobvY6xYkVyemjo1SRodfn8Dql+q+iVARPwS2BMRR9LlAF7OcnJJCyTtkLRL0vIy23slbUlfGyWdm66fImmDpO2Stkn6VMkxp0m6X9IT6fvrc7TXzKx+enth9244ejR574KkAdV7HBMk/Qmgks+ky2eOdmJJY4GbgEuAvcAjktZFxOMluz0FzI2Ig5IuBfqAC4DDwGci4seSxgOPSro/PXY58GBE3JAmo+XAtXkabWZmtavW47gFGA+cUvJ5eHlNhnOfD+yKiCcj4lfAXcDC0h0iYmNEHEwXNwGT0/X7IuLH6ecXgO3ApHS/hcDt6efbgfdniMXMzOqk2nDc6yttk/S7Gc49CXi6ZHkvSW+ikquAe8r8rOkkN+R/lK7qiYh9aYz7JE2oEONSkpIp9PT0MDAwkCHk1jI0NNSWcdeq29oLbnO36LQ2Z54BUNJM4ApgEfBzYM5oh5RZFxXOPZ8kcbx9xPpTgG8Bn46I57PGChARfSSXvpgzZ07Mmzcvz+EtYWBggHaMu1bd1l5wm7tFp7W5auKQNI0kUSwiue8wDZgTEbsznHsvMKVkeTLwbJmfMYvk0telEfFcyfrjSZJG/4jyJvslTUx7GxOBAxliMTOzOqn2HMdGYD1wPPCHETEbeCFj0gB4BJgh6SxJJ5D0VtaN+BlTSUqXLI6InSXrBdwKbI+IL4847zrgyvTzlcB3MsZjZmZ1UO3m+E9Jbob38OooqrKXmsqJiMPAMuBekpvb34yIbZKulnR1utt1JGXbb5b0WDrfB8DbgMXARen6xyRdlm67AbhE0hMkI7ZuyBqTmZkdu2o3xxdK+k3gg8D1kt5AUjH3/Ij4hywnj4j1JL2W0nWrSz4vAZaUOe4HlL9HQno56+IsP9/MLJf+/q58oC+vqvc4IuLnJOVF1krqAT4M/HdJUyJiSrVjzczaSheXSc8r85zjEbE/Ir4aEb/HiNFPZmZtr1qZdHuNajfHPy5pRvpZkv6HpJ9L2gKc1rAIzcwaoYvLpOdVrcfxKWB3+nkRMAs4G/gT4KvFhmVm1mBdXCY9r2qJ43BEDBczfA/wPyPiuYh4gKTsiJlZ5+jiMul5VUscRyVNTMuqXww8ULLtpGLDMjNrsN5e6OuDadNASt77+nxjvIxqo6quAzYDY4F1EbENklkBgScbEJuZWWP19jpRZFDtOY6/TUuOjC+pYAtJMvlw4ZGZmVlLGu05jsPAwRHrXiw0IjMza2mZn+MwMzODURJH+vyGnxA3M7NXVE0c6fzif9OYUMzMrB1kuVS1KeOMf2Zm1gWyzAA4H7ha0m7gRZKqtRERs4oMzMzMWlOWHselJKVGLgLeS/IU+XuLDMrM7Jj198P06TBmTPLe39/siDrGqIkjIvaQTAF7Ufr5pSzHmZk1zXCJ9D17IOLVEulOHnUxagKQ9GfAtcBn01XHA3cUGZSZ2TFxifRCZek5XA68j+T+BhHxLMmUsmZmrckl0guVJXH8Kh2WGwCSTi42JDOzY+QS6YXKkji+KekvSeYb/zhJldxbig3LzOwYuER6oUYdjhsRN0q6BHgeOAe4LiLuLzwyM7NaDVe4XbEiuTw1dWqSNFz5ti6yPMdBmiicLMysfbhEemEqJg5JL5De1ygnIk4tJCIzM2tp1ebjGA8g6c+BfwG+QfLUeC8eVWVm1rWy3Bx/d0TcHBEvRMTzEfF14INFB2ZmZq0pS+I4IqlX0lhJYyT1AkeKDszMzFpTlsTxEeCPgP3p60PpOjMz60JZalXtjoiFEXFGRJwZEe+PiN1ZTi5pgaQdknZJWl5me6+kLelro6RzS7atlXRA0tYRx3xe0jOSHktfl2WJxczM6iNLraozJX1OUl/6Zb5W0toMx40FbiKprjsTWCRp5ojdngLmpiXavwD0lWy7DVhQ4fSrIuK89LV+tFjMrM250m1LyfIcx3eAvyd5YjzPvY3zgV0R8SSApLuAhcDjwztExMaS/TcBk0u2PSxpeo6fZ2adaLjS7XDRwuFKt+DnNJokS+IYFxHX1nDuScDTJct7gQuq7H8VcE/Gcy+T9FFgM/CZiDhYQ3xm1g6qVbp14miKLInjbyVdVsMlIZVZV/aBQknzSRLH2zOc9+skl7Uiff8S8LEy51wKLAXo6elhYGAgU9CtZGhoqC3jrlW3tRfc5izmDg6W/zIZHOShNvlv13G/54io+gJeAI4CvyCpV/UC8HyG4y4E7i1Z/izw2TL7zQL+GXhjmW3Tga1VfkbV7cOv2bNnRzvasGFDs0NoqG5rb4TbnMm0aRHJdEyvfU2bVkB0xWjX3zOwOcp8p2YZVTU+IsZExG9ExKnpcpZyI48AMySdJekE4ApgXekOkqYCdwOLI2JnhnMiaWLJ4uXA1kr7mlkHcKXbljPqpSpJ7yi3PiIernZcRByWtAy4FxgLrI2IbZKuTrevBq4DTgdulgRwOCLmpD/3TmAecIakvcCfRcStwBclnUdyqWo38InRm2lmbcuVbltOlnsc/6nk80kko6UeBS4a7cBI7ousH7FudcnnJcCSCscuqrB+8eghm1lHcaXblpJlPo73li5LmgJ8sbCIzMyspWUpOTLSXuDN9Q7EzMzaQ5Z7HF/j1WG0Y4DzgH8sMCYzM2thWe5xbC75fBi4MyL+T0HxmJlZi8tyj+P2RgRiZmbtoZZ7HGZm1sWcOMys8Vzttq1lucdhZlY3Ex54AFatcrXbNlYxcUj6LhWKEgJExPsKicjMOtrZa9a42m2bq9bjuLFhUZhZ1zjxwIHyGwYHGxuI1axi4oiIh4Y/p0UK35gu7oiIl4sOzMw606EJEzhp//5f3zB1auODsZpkmTp2HvAEyTSwNwM7KxU+NDMbzZNLlrjabZvLMqrqS8C7ImJuRLwDeDewqtiwzKxTHXjnO6GvD6ZNAyl57+vz/Y02kmVU1fERsWN4ISJ2Sjq+wJjMrNO52m1by1RyRNKtwDfS5V6SsupmZtaFKl6qkjQn/fgfgG3AJ4FPAY8DVxcfmpmZtaJqPY5bJJ0C3AncFRFfblBMZmbWwir2OCLiLcB7gCPAX0t6TNK1kqY1LDozM2s5VUdVRcSOiLg+ImYCVwKvA/5Oksuqm5l1qUxFDiWNASYAPcDJwE+LDMrM2ogLFnadqqOqJP0+sAh4P7AVuAu4JiJ+XnxoZtby+vuTAoUuWNhVqo2qehq4AdgOvCUi3hURa500zOwVK1ZULlhoHataj+PtEbGnYZGYWfupVJjQBQs7WrVRVU4aZlZdpcKELljY0TwDoJnVbuVKFyzsQk4cZla73l4XLOxCo9aqknQ28BXgQuAo8EOSkVVPFhybmbUDFyzsOll6HP8L+CbwW8BvA39FUobEzMy6UJbEoYj4RkQcTl93UGUu8tccKC2QtEPSLknLy2zvlbQlfW2UdG7JtrWSDkjaOuKY0yTdL+mJ9P31WWIxM7P6yJI4NkhaLmm6pGmS/jPwvfQL/LRKB0kaSzJr4KXATGCRpJkjdnsKmBsRs4AvAH0l224DFpQ59XLgwYiYATyYLpuZWYNkmY/jw+n7J0as/xhJz+PsCsedD+wavhci6S5gIUlZdgAiYmPJ/puAySXbHpY0vcx5FwLz0s+3AwPAtaM3w8zM6mHUxBERZ9V47knA0yXLe4ELqux/FXBPhvP2RMS+NLZ9kiaU20nSUmApQE9PDwMDA1libilDQ0NtGXetuq290HptnvDAA5y9Zg0nHjjAoQkTeHLJkmSq1zpqtTY3Qqe1OcuoqrHAHwDTS/fPMD+Hyqwre29E0nySxPH20eLJKiL6SC99zZkzJ+bNm1evUzfMwMAA7Rh3rbqtvdBibe7vh1WrXikhctL+/cxctYqZb3pTXUdNtVSbG6TT2pzlHsd3gX8HnA6ML3mNZi8wpWR5MvDsyJ0kzQLWAAsj4rkM590vaWJ67ETgQIZjzGw0rjtlGWW5xzE5vXmd1yPADElnAc8AVwAfKd1B0lTgbmBxROzMeN51JHOD3JC+f6eG2MxsJNedsoyy9DjukfSuvCeOiMPAMuBekgq734yIbZKuljQ8Z/l1JD2Zm9MZBjcPHy/pTpKHDc+RtFfSVemmG4BLJD0BXJIum9mxct0pyyhLj2MT8O10MqeXSe5dREScOtqBEbEeWD9i3eqSz0uAJRWOXVRh/XPAxRniNrM8Vq587dwa4LpTVlaWHseXSMqNjIuIUyNifJakYWZtxnWnLKMsPY4ngK0RkelpcTNrY647ZRlkSRz7gAFJ9wCHhldmGI5rZmYdKEvieCp9nZC+zMysi2V5cvz6RgRiZmbtoWLikPRdqlTBjYj3FRKRmZm1tGo9jhsbFoWZFaO/P3nye3AweR5j5Urf/LZjVjFxRMRDjQzEzOqsv/+1z2Xs2ZMsg5OHHRPPOW7WqVx7ygrixGHWqVx7ygrixGHWqVx7ygpSbVTVumoHelSVWYtz7SkrSLVRVReSzOB3J/Ajyk/MZGatavgGuEdVWZ1VSxy/RVK2fBHJPBrfA+6MiG2NCMzM6sC1p6wAFe9xRMSRiPh+RFwJvBXYRVKz6o8bFp2ZmbWcqiVHJJ1IMt/4IpI5x79KMmOfmZl1qWo3x28H3gzcA1wfEVsbFpWZmbWsaj2OxcCLwBuBT0qv3BvPPAOgmZl1nmr3OMaks/2NT2f+O9UzAJo1UX8/TJ8OY8Yk7/39zY7IulS1S1XjgJcj4uV0+RzgMmB3RHy7QfGZGbjulLWUak+Of5/khjiS3gD8EDgbWCbphuJDM7NXuO6UtZBqieP1EfFE+vlKkmc4/hi4lGSklZk1iutOWQupljhKJ3G6CLgfICJ+BRwtMigzG8F1p6yFVEscWyTdKOka4A3AfQCSXteIwMysxMqVSZ2pUq47ZU1SLXF8HPgZyX2Od0XE8AXWmXh2QLPG6u2Fvj6YNg2k5L2vzzfGrSmqzQD4C+A1N8ElnQH8MCI2Fh2YmY3gulPWIir2OCS9VdKApLslvUXSVmArsF/SgsaFaGZmraTak+N/AXwO+E3g74BLI2KTpN8hKbX+/QbEZ2ZmLabaPY7jIuK+iPgr4F8iYhNARPwk68klLZC0Q9IuScvLbO+VtCV9bZR07mjHSvq8pGckPZa+Lssaj1nLSJ8Cn3vRRX4K3NpOtR5H6ZDbX4zYFoxC0ljgJpI5PfYCj0haFxGPl+z2FDA3Ig5KuhToAy7IcOyqiPANemtPJU+BC/wUuLWdaj2OcyU9L+kFYFb6eXj5X2c49/nAroh4Mn324y5gYekOEbExIg6mi5uAyVmPNWtbfgrc2ly1UVVjj/Hck0imnh22F7igyv5XkZRwz3LsMkkfBTYDnylJPq+QtBRYCtDT08PAwEDe+JtuaGioLeOuVbe0d+7gYNl5mGNwkIe6oP3d8nsu1WltrjqR0zEq+/9G2R2l+SSJ4+0Zjv068IV0+QvAl4CP/drOEX0kl76YM2dOzJs3L0forWFgYIB2jLtWXdPeqVOTy1MjaOrUrmh/1/yeS3Ram6tdqjpWe4EpJcuTgWdH7iRpFrAGWBgRz412bETsT6e1PQrcQnJZy6x9+Clwa3NFJo5HgBmSzpJ0AnAFsK50B0lTSaaiXRwRO7McK2liyX6XkzxbYtY+Sp4CDz8Fbm2osEtVEXFY0jLgXmAssDYitkm6Ot2+GrgOOB24OZ1h8HBEzKl0bHrqL0o6j+RS1W7gE0W1waww6VPgD3XYJQzrDkXe4yAi1gPrR6xbXfJ5CbAk67Hp+sV1DrO8/v5klMvgYHJNeuVK/0VoZkaxl6ra1/A4+z17IOLVcfZ+SMsq8bSu1kWcOMrxOHvLw39oWJdx4ijHs61ZHv5Dw7qME0c5nm3N8vAfGtZlnDjK8Th7y8N/aFiXceIox7OtWR7+Q8O6TKHDcduaZ1uzrIb/nXj4tnUJJw6zevAfGtZFfKmqnjyWvzP492hWlXsc9VIyOQ/gyXnalX+PZqNyj6NePJa/M/j3aDYqJ4568Vj+zuDfo9monDjqxWP5O4N/j2ajcuKoF4/l7wz+PZqNyomjXvzQYGfw79FsVE4c9dTbC7t3w9Gjybu/bJqvlqG1/j2aVeXE0Wx+ZqA4LnduVggnjmbyF1uxPLTWrBBOHM3kL7ZieWitWSGcOJrJX2zF8tBas0I4cTSTv9iK5aG1ZoVw4mgmf7Hll2cwgYfWmhXCiaOZavli6+ZRWLUMJvDQWrO6c+JotjxfbN0+CsuDCcxaghNHO+m0L84RvacJDzxQfX8PJjBrCU4c7aTWL85WvLxVpvd0zo03Vo/NgwnMWoITRzup5YuzVS9vlek9jT10qHrvyYMJzFqCE0c7qeWLs9bLW3l7KXn3r6X35FFSZi2h0MQhaYGkHZJ2SVpeZnuvpC3pa6Okc0c7VtJpku6X9ET6/voi29BSavnirOULOm8vpZZeTa2XnTxKyqzpCkscksYCNwGXAjOBRZJmjtjtKWBuRMwCvgD0ZTh2OfBgRMwAHkyXu0feL85avqDz9lJq6dWU6T0dOfFEX3YyawNF9jjOB3ZFxJMR8SvgLmBh6Q4RsTEiDqaLm4DJGY5dCNyefr4deH9xTegAtVzeyttLqdNlpx1/+qfuQZi1geMKPPck4OmS5b3ABVX2vwq4J8OxPRGxDyAi9kmaUO5kkpYCSwF6enoYGBjIG3/TDQ0NHXvckyYx4ZprOHvNGk48cIBDEybw5JIlHJg0CSqc+60TJnDS/v2/tv6XEyawqcwxefcvjY3bbntlcWhoiANt+Hs6FnX5HbcZt7kDREQhL+BDwJqS5cXA1yrsOx/YDpw+2rHA/xtx7MHRYpk9e3a0ow0bNjTnB99xR8S4cRHJHYvkNW5csr4e+1fQtPY2kdvcHdq1zcDmKPOdWuSlqr3AlJLlycCzI3eSNAtYAyyMiOcyHLtf0sT02InAgTrHbXlvwnu0k1lXKfJS1SPADElnAc8AVwAfKd1B0lTgbmBxROzMeOw64ErghvT9OwW2oXv19ub74s+7v5m1rcISR0QclrQMuBcYC6yNiG2Srk63rwauA04HbpYEcDgi5lQ6Nj31DcA3JV0FDJJc1jIzswYpssdBRKwH1o9Yt7rk8xJgSdZj0/XPARfXN1IzM8vKT46bmVkuThxmZpaLE4eZmeWiZKhuZ5P0U2BPs+OowRnAz5odRAN1W3vBbe4W7drmaRFx5siVXZE42pWkzRExp9lxNEq3tRfc5m7RaW32pSozM8vFicPMzHJx4mhtfc0OoMG6rb3gNneLjmqz73GYmVku7nGYmVkuThxmZpaLE0eLkHSNpG2Stkq6U9JJnT6/uqRPpe3dJunT6bqOarOktZIOSNpasq5iGyV9VtIuSTskvbs5UR+bCm3+UPp7Pippzoj9O7XN/03STyRtkfRtSa8r2dbWbXbiaAGSJgGfBOZExJtJKgJfQQfPry7pzcDHSaYJPhd4j6QZdF6bbwMWjFhXto2SZpL83v9VeszNksY2LtS6uY1fb/NW4APAw6UrO7zN9wNvjohZwE7gs9AZbXbiaB3HAb8h6ThgHMnEVZ08v/qbgE0R8VJEHAYeAi6nw9ocEQ8D/3fE6kptXAjcFRGHIuIpYBdJYm0r5docEdsjYkeZ3Tu5zfel/7YBNpFMSAcd0GYnjhYQEc8AN5LML7IP+HlE3MeI+dWBsvOrt6mtwDsknS5pHHAZyayPndzmYZXaOAl4umS/vem6TtYtbf4YcE/6ue3b7MTRAtJr3AuBs4DfBk6W9G+bG1WxImI78F9JuvPfB/4ROFz1oM6nMus6fbx8x7dZ0gqSf9v9w6vK7NZWbXbiaA3vBJ6KiJ9GxMsk0+n+Hh0+v3pE3BoR/yYi3kHSzX+CDm9zqlIb95L0uoZNJrlk2ck6us2SrgTeA/TGqw/NtX2bnThawyDwVknjlMyhezGwnVfnV4cOnF9d0oT0fSrJjdM76fA2pyq1cR1whaQTJZ0FzAD+oQnxNVLHtlnSAuBa4H0R8VLJpvZvc0T41QIv4HrgJyTX/r8BnEgyH/uDJH+JPwic1uw469zmvwceJ7lMdXG6rqPaTJIM9wEvk/yleVW1NgIrgH8GdgCXNjv+Orb58vTzIWA/cG8XtHkXyb2Mx9LX6k5ps0uOmJlZLr5UZWZmuThxmJlZLk4cZmaWixOHmZnl4sRhZma5OHGY1ZGkFWkV2C2SHpN0gaRPp2VVKh2zJi18h6ShxkVrVhsPxzWrE0kXAl8G5kXEIUlnACcAG0kqH/+szDFjI+JIyfJQRJzSsKDNauAeh1n9TAR+FhGHANJE8Yck9cc2SNoASXKQ9OeSfgRcKGmgzBwVZ0j6oaQ/kHSmpG9JeiR9va3B7TJ7DScOs/q5D5giaaekmyXNjYivktQhmh8R89P9Tga2RsQFEfGDkSeR1AN8D7guIr4HfAVYFRG/C3wQWNOQ1phVcFyzAzDrFBExJGk28PvAfOB/Syo3EdUR4FsVTnM8SRmS/xgRD6Xr3gnMTMqYAXCqpPER8UL9ojfLzonDrI7S+xUDwICkf+LVYoalfll6X2OEw8CjwLtJJreC5MrAhRHxizqHa1YTX6oyqxNJ56TT3w47D9gDvACMz3iaIJn053dKeiv3ActKfs55xxys2TFwj8Osfk4BvibpdSQ9h13AUmARcI+kfSX3OSqKiCOSrgC+K+l5kvnob5K0heT/2YeBqwtqg9moPBzXzMxy8aUqMzPLxYnDzMxyceIwM7NcnDjMzCwXJw4zM8vFicPMzHJx4jAzs1z+Pw2RSYedcylaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# volatility smile plot: try with different parameters\n",
    "plt.plot(strike, bsm_vol, 'ro')\n",
    "#plt.axis([0, 6, 0, 20])\n",
    "plt.ylabel('BSM Imp Vol under SABR')\n",
    "plt.xlabel('Strike')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (2) SABR for $\\beta=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "strike = np.linspace(75,125,num=25)\n",
    "#strike = 1\n",
    "forward = 100\n",
    "sigma = 20\n",
    "texp = 1\n",
    "alpha = 0.5\n",
    "rho = -0.25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sabr_norm = pf.SabrHagan2002(sigma, vov=vov, rho=rho, beta=0)\n",
    "price = sabr_norm.price(strike, forward, texp)\n",
    "nvol = sabr_norm.vol_smile(strike, forward, texp, model='norm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEJCAYAAACOr7BbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfmUlEQVR4nO3dfZRcdZ3n8feH8GTScRUwPRDyAMogGQaEtIjDrJNGnswyoIzsCRvZLKI57gLiCDOA8fgwMznLjgHXx9UeYAalsdUDrLhESMA0DI4gCUZIbAIRSIyJYfCJtEgwyXf/uLel0rlVfavStx4/r3PqVN1f3Vv1/aWT/ub+HhURmJmZjbZPowMwM7Pm5ARhZmaZnCDMzCyTE4SZmWVygjAzs0xOEGZmlqmwBCFpmqQVkoYkrZV0eVr+95Iek7Ra0jJJh5W5/ixJ6yStl3R1UXGamVk2FTUPQtKhwKER8aikycAq4J3Apoh4IT3ng8CsiPjAqGsnAE8CpwObgEeACyLix4UEa2Zme9i3qA+OiC3AlvT1NklDwNRRv+QnAVkZ6iRgfUQ8DSBpADgXqJggDjnkkJg5c+Y4RF8/v/3tb5k0aVKjw6gr17kzuM6tYdWqVc9HxOuy3issQZSSNBM4AXg4PV4M/FfgN0BvxiVTgZ+WHG8C3jLW98ycOZOVK1fubbh1NTg4yJw5cxodRl25zp3BdW4NkjaUfa/opTYkdQH3A4sj4vZR710DHBgRHx9Vfj5wZkS8Lz2+EDgpIi7L+PyFwEKA7u7u2QMDA8VUpCDDw8N0dXU1Ooy6cp07g+vcGnp7e1dFRE/We4XeQUjaD7gN6B+dHFK3AncBHx9VvgmYVnJ8OLA56zsiog/oA+jp6YlWy96t+D+OveU6dwbXufUVOYpJwI3AUERcX1J+VMlp5wBPZFz+CHCUpCMk7Q/MA+4sKlYzM9tTkXcQpwAXAo9LWp2WfQS4WNLRwC5gA/ABgHS46w0RMTcidki6FLgHmADcFBFrC4zVzMxGKXIU04OAMt5aWub8zcDckuOl5c41M7PieSZ1fz/MnAn77JM89/c3OiIzs6ZQl2GuTau/HxYuhBdfTI43bEiOAebPb1xcZmZNoLPvIBYteiU5jHjxxaTczKzDdXaC2LixunIzsw7S2Qli+vTqys3MOkhnJ4jFi2HixN3LJk5Mys3MOlxnJ4j586GvD2bMACl57utzB7WZGZ0+igmSZOCEYGa2h86+gzAzs7KcIMzMLJMThJmZZXKCMDOzTE4QZmaWyQnCzMwyOUGYmVkmJwgzM8vkBGFmZpmcIKrlDYbMrEN4qY1qeIMhM+sgvoOohjcYMrMO4gRRDW8wZGYdxAmiGt5gyMw6iBNENbzBkJl1ECeIaniDITPrIIWNYpI0DfgK8EfALqAvIj4j6VPAXwIvAz8BLoqIX2dc/yywDdgJ7IiInqJirYo3GDKzDlHkHcQO4IqIOAY4GbhE0ixgOXBsRBwHPAlcU+EzeiPiTU2THMzMOkhhCSIitkTEo+nrbcAQMDUilkXEjvS0h4DDi4rBzMxqV5c+CEkzgROAh0e99V7gO2UuC2CZpFWSFhYYnpmZZVBEFPsFUhdwP7A4Im4vKV8E9ADnRUYQkg6LiM2SppA0S10WEQ9knLcQWAjQ3d09e2BgoKCaFGN4eJiurq5Gh1FXrnNncJ1bQ29v76pyzfiFJghJ+wH/D7gnIq4vKV8AfAB4e0S8WO76kvM/AQxHxJJK5/X09MTKlSv3Lug6GxwcZM6cOY0Oo65c587gOrcGSWUTRGFNTJIE3AgMjUoOZwFXAeeUSw6SJkmaPPIaOANYU1SsZma2pyL7IE4BLgROlbQ6fcwFPg9MBpanZV+CpElJ0tL02m7gQUk/An4A3BURdxcYq5mZjVLYPIiIeBBQxltLM8qIiM3A3PT108DxRcVmZmZj80xqMzPLVFOCkHTKeAfS1rzJkJm1oLJNTJImAP8ZmArcHRFrJJ0NfAR4Fcm8BhuLNxkysxZV6Q7iRuB9wMHAZyX9M7AE+MeIcHLIy5sMmVmLqtRJ3QMcFxG7JB0IPA+8ISJ+Xp/Q2oQ3GTKzFlXpDuLliNgFEBEvAU86OdTAmwyZWYuqlCDeKOmx9PF4yfHjkh6rV4Atz5sMmVmLqtTEdEzdomhnIx3RixYlzUrTpyfJwR3UZtbkyiaIiNgwukzSIcAvshbXswq8yZCZtaCyTUySTpY0KOl2SSdIWkOyHtLWdD0lMzNrY5WamD5PMufhPwDfBd4REQ9JeiPwNcBrI5mZtbFKndT7pru/fRP4eUQ8BBART9QnNDMza6RKCWJXyevfjXrPfRBmZm2uUhPT8ZJeIFmR9VXpa9LjAwuPzMzMGqrSKKYJ9QzEzMyaS1WruaY7vc2XdFdRAZmZWXMYM0FI2l/SOyV9A9gCnAZ8qfDIzMysoSrNgzhd0k3AM8C7ga8Cv4yIiyLi2/UKsGN5Dwkza7BKndT3AP8K/HlEPAMg6TN1iarTeQ8JM2sClZqYZgMPAfdKWi7pYsAd1/XgPSTMrAmUTRAR8cOIuCoiXg98gmQHuf0lfUfSwnoF2JG8h4SZNYFco5gi4nsRcSnJ9qOfBt5aaFSdzntImFkTqNRJPVHSfiXHRwOXA10RcVE9gutY3kPCzJpApTuIu4GZAJLeAHwfOBK4RNL/LD60DjZ/PvT1wYwZICXPfX3uoDazuqo0ium1EfFU+noB8LWIuEzS/sAq4JrCo+tk3kPCzBqs0h1E6YJ8pwLLASLiZXZfyC+TpGmSVkgakrRW0uVp+ackPZFuX3qHpNeUuf4sSeskrZd0de4amZnZuKiUIB6TtETSXwNvAJYBlPuFnmEHcEVEHAOcTNI0NYsk0RwbEccBT5JxJyJpAvAF4B3ALOCC9FozM6uTSgni/cDzJP0QZ0TEyMD8WcCSsT44IrZExKPp623AEDA13WNiR3raQ8DhGZefBKyPiKfTO5YB4Nwc9TEzs3GiemwvLWkm8ADJncMLJeXfBr4eEbeMOv/dwFkR8b70+ELgLelQ29GfvRBYCNDd3T17YGCgsHoUYXh4mK6urkaHUVeuc2dwnVtDb2/vqojoyXqvUif1uJDUBdwGfGhUclhE0gyVtciQMsoyM1lE9AF9AD09PTFnzpy9DbmuBgcHabWY95br3Blc59ZXaIJI51HcBvRHxO0l5QuAs4G3R/YtzCZgWsnx4cDmImM1M7PdVZxJLWmCpE/V8sGSBNwIDEXE9SXlZwFXAeeU9GuM9ghwlKQj0mG184A7a4nDzMxqUzFBRMROYHb6y75apwAXAqdKWp0+5gKfByYDy9OyLwFIOkzS0vR7dwCXkqwoOwR8IyLW1hCDmZnVKE8T0w+Bb0n6JvDbkcLSJqMsEfEg2X0JS8ucvxmYW3K8tNy5VkZ/f7Li68aNybpNixd7sp2Z1SxPgjgI+AXJZLkRAVRMEFZn3kPCzMbZmAnCC/O1iEp7SDhBmFkN8uxJ/ceS7pO0Jj0+TtJHiw/NquI9JMxsnOXZD+KfSJbD+D1ARDxGMqrImon3kDCzcZYnQUyMiB+MKtuReaY1jveQMLNxlidBPC/p9aQzmdNlMLYUGpVVz3tImNk4yzOK6RKSpSzeKOlnwDPAewqNymrjPSTMbBzlGcX0NHCapEnAPunKrGZm1ubKJghJHy5TDkDp8hlmZtZ+Kt1BTE6fjwbezCtrIf0lydLdZmbWxsomiIj4JICkZcCJI01Lkj4BfLMu0ZmZWcPkGcU0HXi55Phlkl3mzMysjeUZxfRV4AeS7iAZ6vou4CuFRmVmZg2XZxTTYkl3A3+eFl0UET8sNiwzM2u0PE1MAKtJ+h3uAH4hyes3tIP+fpg5E/bZJ3nuz9r91cw61Zh3EJIuAz4ObAV2kuzxEMBxxYZmhfLy4GY2hjx3EJcDR0fEn0TEcRHxpxHh5NDqKi0PbmZGvgTxU+A3RQdideblwc1sDHlGMT0NDEq6C9g+UuiZ1C1u+vSkWSmr3MyMfHcQG4HlwP4ks6tHHtbKvDy4mY0hzzDXT9YjEKuzkY7oRYuSZqXp05Pk4A5qM0vlGcW0gnQviFIRcWohEVn9eHlwM6sgTx/ElSWvDwT+Cu8oZ2bW9vI0Ma0aVfQ9SfcXFI+ZmTWJMTupJR1U8jhE0pnAH+W4bpqkFZKGJK2VdHlafn56vEtST4Xrn5X0uKTVklZWVSszM9treZqYVpH0QYikaekZ4OIc1+0AroiIRyVNBlZJWg6sAc4DvpzjM3oj4vkc55mZ2TjL08R0RC0fHBFbgC3p622ShoCpEbEcXtmZzszMmpMi9higNP5fIs0k2YXu2Ih4IS0bBK6MiMzmI0nPAL8iuXv5ckT0lTlvIbAQoLu7e/bAwMC4x1+k4eFhurq6Gh1GXbnOncF1bg29vb2rIiKzuT9PE9NekdQF3AZ8aCQ55HRKRGyWNAVYLumJiNhjq9M0cfQB9PT0xJw5c8Yj7LoZHByk1WLeW65zZ3CdW1/e5b5rImk/kuTQHxG3V3NtRGxOn58jWWb8pPGP0KrmJcLNOkbZOwhJJ1a6MCIerfS+kk6GG4GhatdtkjQJ2Cftu5gEnAH8XTWfYQXwEuFmHaVSE9N1Fd4LYKyZ1KcAFwKPS1qdln0EOAD4HPA64C5JqyPiTEmHATdExFygG7gj7cjeF7g1Iu4eqzJWsEpLhDtBmLWdsgkiInr35oMj4kGSobFZ7sg4fzMwN339NHD83ny/FcBLhJt1lDxrMe0H/HfgbWnRIMmoot8XGJc1Iy8RbtZR8nRS/x9gNvDF9DE7LbNO4yXCzTpKnmGub46I0uae70r6UVEBWRPzEuFmHSVPgtgp6fUR8RMASUcCO4sNy5qWlwg36xh5EsTfACskPU3S6TwDuKjQqMzMrOEqzYO4C7gV+L/AUcDRJAniiYjYXu46MzNrD5U6qfuAs0lWb/0KcCTJpDcnBzOzDlA2QUTEtyLiApImpduBBcBGSTdJOr1eAZqZWWOMOcw1In4XEV+PiHeRLHlxAuBZzWZmbS7PjnLdki6T9D2S/ohlJHMhzMysjZVNEJLeL+m7wKPAHwN/GxFHRsRVEbG6XgFaGxi1AuyUe+9tdERmlkOlYa5/BlwL3BsRu+oUj7WbjBVgj16yBI45xvMpzJpcpU7qiyJimZOD7ZWMFWAnbN+elJtZUyt0wyAzrwBr1rqcIKxY5VZ69QqwZk0vV4KQdKKkD6ajmSruNGe2m4wVYHcecIBXgDVrAXmGuX4MuBk4GDgE+GdJHy06MGsT8+dDXx/MmAESzJjBuiuvdAe1WQvIs1jfBcAJEfESgKRrSYa+/kORgVkbGbUC7HODg8xqYDhmlk+eJqZngQNLjg8AflJINGZm1jTy3EFsB9ZKWg4EcDrwoKTPAkTEBwuMz8zMGiRPgrgjfYwYLCYUMzNrJmMmiIi4uR6BmJlZc8kziulsST+U9EtJL0jaJumFegRnZmaNk6eT+n+T7AVxcES8OiImR8Sriw3LOtqoxf3o7290RGYdKU+C+CmwJiKimg+WNE3SCklDktZKujwtPz893iWpp8L1Z0laJ2m9pKur+W5rYSOL+23YABHJ88KFThJmDZCnk/pvgaWS7icZ0QRARFw/xnU7gCsi4lFJk4FV6UioNcB5wJfLXShpAvAFkhFTm4BHJN0ZET/OEa+1sozF/XjxxaTck+vM6ipPglgMDJPMhdg/7wdHxBZgS/p6m6QhYGpELAeQVOnyk4D1EfF0eu4AcC7gBNHuvLifWdPIkyAOiogz9uZLJM0k2ar04ZyXTCVp2hqxCXhLmc9eCCwE6O7uZnBwsOY4G2F4eLjlYt5blep88pQpHLh16x7lL02ZwkMt/Ofkn3NnaLc650kQ90o6IyKW1fIFkrqA24APRUTe0U9ZtxeZfSAR0Qf0AfT09MScOXNqCbNhBgcHabWY91bFOl933e4bDAFMnMiB113X0n9O/jl3hnarc55O6kuAuyX9rtphrpL2I0kO/RFxexVxbQKmlRwfDmyu4nprVRmL+9HX5/4HswbIM1Fuci0frKST4UZgKEeH9miPAEdJOgL4GTAP+C+1xGEtaNTifmbWGGUTxFj7PkTEo2N89inAhcDjklanZR8hWezvc8DrgLskrY6IMyUdBtwQEXMjYoekS4F7gAnATRGxNleNzMxsXFS6g7iuwnsBnFrpgyPiQbL7EmD3tZ1Gzt8MzC05XgosrfQdZmZWnLIJIiJ66xmImZk1F+9JbWZmmZwgrD14/SazcZdnHoRZcxtZv2lk7sTI+k3g0VBme6HIUUxm9eH1m8wKUdgoJrO68fpNZoXwKCZrfdOnJ81KWeVmVrNcfRCSjgVmkazoCkBEfKWooMyqsnhx5vpNLF7cuJjM2kCeLUc/TjLz+XNAL/CPwDkFx2WWn9dvMitEnjuIdwPHAz+MiIskdQM3FBuWWZW8fpPZuMszD+J3EbEL2CHp1cBzwJHFhmVmZo2W5w5ipaTXAP8ErCLZXe4HRQZlZmaNl2e57/+RvvySpLuBV0fEY8WGZWZmjZZrqQ1Jx0k6BzgReIOk84oNy6wOvDyHWUVj3kFIugk4DlgL7EqLA6hmhziz5uLlOczGlKcP4uSImFV4JGb15OU5zMaUp4np+5KcIKy9eHkOszHluYO4mSRJ/BzYTrJLXETEcYVGZlYkL89hNqY8CeIm0r2leaUPwqy1eXkOszHlSRAbI+LOwiMxq6eRfoZFi5JmpenTk+Tg/gezP8iTIJ6QdCvwbZImJgAiwqOYrLV5eQ6zivIkiFeRJIYzSso8zNXMrM1VTBCSJgDPR8Tf1CkeMzNrEhWHuUbETpLZ02Zm1mHyzINYLelOSRdKOm/kMdZFkqZJWiFpSNJaSZen5QdJWi7pqfT5tWWuf1bS45JWS1pZZb3Mxp+X5rAOk6cP4iDgF+y+B3WePogdwBUR8aikycAqScuB/wbcFxHXSroauBq4qsxn9EbE8zliNCuWl+awDpRnNdeLavngiNgCbElfb5M0BEwFzgXmpKfdDAxSPkGYNQcvzWEdKM+Wo4dLukPSc5K2SrpN0uHVfImkmcAJwMNAd5o8RpLIlDKXBbBM0ipJC6v5PrNx56U5rAMpIiqfkDQL3Qp8NS16DzA/Ik7P9QVSF3A/sDgibpf064h4Tcn7v4qIPfohJB0WEZslTQGWA5dFxAMZ5y0EFgJ0d3fPHhgYyBNW0xgeHqarq6vRYdRVK9b55HnzOHDr1j3KX+ru5qEcf+dasc57y3VuDb29vasioifzzYio+ABW5ykrc+1+wD3Ah0vK1gGHpq8PBdbl+JxPAFeOdd7s2bOj1axYsaLRIdRdS9b5llsiJk6MgFceEycm5Tm0ZJ33kuvcGoCVUeZ3ap5RTM9Leo+kCenjPSSd1hVJEnAjMBQR15e8dSewIH29APhWxrWT0o5tJE0imaS3JkesZsWYPx/6+mDGDJCS574+9z9YW8szium9wOeBT5P0C/xbWjaWU0gX+ZO0Oi37CHAt8A1JFwMbgfMhaVICboiIuUA3cEeSY9gXuDUi7s5ZJ7NieGkO6zB5RjFtBM6p9oMj4kGSpcGzvD3j/M3A3PT108Dx1X6nmZmNn7IJQtLHKlwXEfH3BcRj1l76+2HRIv7CK8ZaC6p0B/HbjLJJwMXAwYAThFklJZPrBJ5cZy2nbCd1RFw38gD6SFZ1vQgYAI6sU3xmravS5DqzFjDWaq4HAR8G5pPMej4xIn5Vj8DMWp4n11mLK3sHIelTwCPANuBPI+ITTg5mVSi3v7X3vbYWUWkexBXAYcBHgc2SXkgf2yS9UJ/wzFrY4sXJPtelvO+1tZCyTUwRkWcSnZmVU7LvdWzciDyKyVqMk4BZkebPh2ef5f7vfheefdbJwVqKE4RZs/HGRNYk8iy1YWb14o2JrIn4DsKsmXjuhDURJwizZuK5E9ZEnCDMmonnTlgTcYIwayaeO2FNxAnCrJl4YyJrIk4QZs0mnTvBrl355k54WKwVxMNczVqZh8VagXwHYdbKPCzWCuQEYdbKPCy2sxXcvOgEYdbKPCy2c400L27YABGvNC+OY5JwgjBrZR4W27nq0LzoBGHWymodFuuRT62vDs2LHsVk1urmz69uxJJHPrWH6dOTn11W+TjxHYRZp/HIp+ZU7V1dHZoXC0sQkqZJWiFpSNJaSZen5QdJWi7pqfT5tWWuP0vSOknrJV1dVJxmHccjn5pPLR3OdZh1X+QdxA7giog4BjgZuETSLOBq4L6IOAq4Lz3ejaQJwBeAdwCzgAvSa81sb3nkU/Op9a6u2ln3VSosQUTEloh4NH29DRgCpgLnAjenp90MvDPj8pOA9RHxdES8DAyk15nZ3qq1acId28Vp0ru6uvRBSJoJnAA8DHRHxBZIkggwJeOSqcBPS443pWVmtrdqaZqow5j7tlNNQm3SuzpFRLFfIHUB9wOLI+J2Sb+OiNeUvP+riHjtqGvOB86MiPelxxcCJ0XEZRmfvxBYCNDd3T17YGCguMoUYHh4mK6urkaHUVeuc+s5ed48Dty6dY/yl7q7eajMv7lWr3MtRuo85d57OXrJEiZs3/6H93YecADrrryS5047bY/rqj1/PPX29q6KiJ7MNyOisAewH3AP8OGSsnXAoenrQ4F1Gde9Fbin5Pga4Jqxvm/27NnRalasWNHoEOrOdW5BUkRy77D7Qyp7ydpFiyJmzEjOmTEj4pZb6hZuo/zh5zxjRvaf14wZ5S++5ZaG/HkBK6PM79QiRzEJuBEYiojrS966E1iQvl4AfCvj8keAoyQdIWl/YF56nZk1QrVNIP39HL1kSfs0SVXb/1JLn0LBHc61KLIP4hTgQuBUSavTx1zgWuB0SU8Bp6fHSDpM0lKAiNgBXEpy9zEEfCMi1hYYq5lVUm3H9qJFuzWXAPlG5TRjR3gt/S9N2qdQtXK3Fq34cBNTa3CdW1Q1TSA1NEnFLbdETJy4+/kTJ47d1FJt00y151fRXPSHn3OtdWkAKjQxNfyX+ng+nCBag+vcAWppg6+13b6aX8S1/OKuItnt9nNuUJ9CtSolCC+1YWbjb/Fidh5wwO5lY821qKXdvtoJZrVMSKu1uagJ+xSq5QRhZuNv/nzWXXlldXMtavlFXG1SqSUJdfCS6k4QZlaI5047rbr/Qdfyi7japFJLEqrDmkfNygnCzJpDLb+Iq00qtd4NtEFzUS2cIMyseVT7i7japNLBdwO18IZBZtbaqt0wqdrzO5jvIMzMLJMThJmZZXKCMDOzTE4QZmaWyQnCzMwyFb5hUD1J+ndgQ6PjqNIhwPONDqLOXOfO4Dq3hhkR8bqsN9oqQbQiSSuj3G5Obcp17gyuc+tzE5OZmWVygjAzs0xOEI3X1+gAGsB17gyuc4tzH4SZmWXyHYSZmWVygqgjSX8taa2kNZK+JulASQdJWi7pqfT5tY2OczxJujyt71pJH0rL2q7Okm6S9JykNSVlZesp6RpJ6yWtk3RmY6KuXZn6np/+nHdJ6hl1fkvXF8rW+VOSnpD0mKQ7JL2m5L2Wr7MTRJ1Imgp8EOiJiGOBCcA84Grgvog4CrgvPW4Lko4F3g+cBBwPnC3pKNqzzv8CnDWqLLOekmaR/Oz/JL3mi5Im1C/UcfEv7FnfNcB5wAOlhW1SX8iu83Lg2Ig4DngSuAbap85OEPW1L/AqSfsCE4HNwLnAzen7NwPvbExohTgGeCgiXoyIHcD9wLtowzpHxAPAL0cVl6vnucBARGyPiGeA9SRJtGVk1TcihiJiXcbpLV9fKFvnZenfbYCHgMPT121RZyeIOomInwFLgI3AFuA3EbEM6I6ILek5W4ApjYty3K0B3ibpYEkTgbnANNq7zqXK1XMq8NOS8zalZe2qU+r7XuA76eu2qLMTRJ2k7c/nAkcAhwGTJL2nsVEVKyKGgP9Fcht+N/AjYEfFizqDMsraeThh29dX0iKSv9v9I0UZp7VcnZ0g6uc04JmI+PeI+D1wO/BnwFZJhwKkz881MMZxFxE3RsSJEfE2ktvzp2jzOpcoV89NJHdSIw4naW5sV21dX0kLgLOB+fHKvIG2qLMTRP1sBE6WNFGSgLcDQ8CdwIL0nAXAtxoUXyEkTUmfp5N0YH6NNq9ziXL1vBOYJ+kASUcARwE/aEB89dK29ZV0FnAVcE5EvFjyVnvUOSL8qNMD+CTwBEnb/FeBA4CDSUa4PJU+H9ToOMe5zv8K/JikeentaVnb1Zkk8W0Bfk/yv8eLK9UTWAT8BFgHvKPR8Y9Tfd+Vvt4ObAXuaZf6VqjzepK+htXp40vtVGfPpDYzs0xuYjIzs0xOEGZmlskJwszMMjlBmJlZJicIMzPL5ARhViVJi9JVSx+TtFrSWyR9KF1OpNw1N6QLuCFpuH7RmtXOw1zNqiDprcD1wJyI2C7pEGB/4N9IVup9PuOaCRGxs+R4OCK66ha0WY18B2FWnUOB5yNiO0CaEN5Nsr7WCkkrIEkCkv5O0sPAWyUNZuyRcIik70v6T5JeJ+k2SY+kj1PqXC+zPThBmFVnGTBN0pOSvijpLyLisyTr7PRGRG963iRgTUS8JSIeHP0hkrqBu4CPRcRdwGeAT0fEm4G/Am6oS23MKti30QGYtZKIGJY0G/iPQC/wdUlZGx7tBG4r8zH7kSy9cUlE3J+WnQbMSpbpAuDVkiZHxLbxi96sOk4QZlVK+xMGgUFJj/PKgnylXirtdxhlB7AKOJNkEyVI7ubfGhG/G+dwzWrmJiazKkg6Ot02dcSbgA3ANmByzo8Jks1l3lhy97EMuLTke96018Ga7SXfQZhVpwv4XLo5/Q6S1TwXAhcA35G0paQfoqyI2ClpHvBtSS+Q7Ff+BUmPkfy7fAD4QEF1MMvFw1zNzCyTm5jMzCyTE4SZmWVygjAzs0xOEGZmlskJwszMMjlBmJlZJicIMzPL5ARhZmaZ/j/+mBvFj+q/agAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(strike, nvol, 'ro')\n",
    "#plt.axis([0, 6, 0, 20])\n",
    "plt.ylabel('Normal Imp Vol under SABR')\n",
    "plt.xlabel('Strike')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3) Smile calibration to 3 options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([13.608,  8.165,  4.769]), array([0.201, 0.205, 0.213]))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "strike3 = np.array([90, 100, 110])\n",
    "price3 = sabr_bsm.price(strike3, forward, texp)\n",
    "vol3 = sabr_bsm.vol_smile(strike3, forward, texp)\n",
    "\n",
    "price3, vol3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'sigma': 0.2,\n",
       " 'intr': 0.0,\n",
       " 'divr': 0.0,\n",
       " 'is_fwd': False,\n",
       " 'vov': 0.5,\n",
       " 'rho': 0.25,\n",
       " 'beta': 1}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sabr_bsm.__dict__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'sigma': 0.2000000000001008,\n",
       " 'vov': 0.4999999998672808,\n",
       " 'rho': 0.24999999998258432}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# makre sure this return the original sigma, alpha, rho\n",
    "sabr_bsm.calibrate3(vol3, strike3, forward, texp, is_vol=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'sigma': 0.19999999979636335,\n",
       " 'vov': 0.5000000129299839,\n",
       " 'rho': 0.2499999894054446}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# makre sure this return the original sigma, alpha, rho\n",
    "sabr_bsm.calibrate3(price3, strike3, forward, texp, is_vol=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Pricing under MC method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "strike = np.linspace(75,125,num=25)\n",
    "forward = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'sigma': 0.2, 'intr': 0.0, 'divr': 0.0, 'is_fwd': False, 'vov': 0.5, 'rho': 0.25, 'beta': 1}\n"
     ]
    }
   ],
   "source": [
    "# instantiate mc model from the hagan model's parameters\n",
    "print(sabr_bsm.__dict__)\n",
    "sabr_bsm_mc = opt.sabr.ModelBsmMC(sabr_bsm.sigma, vov=sabr_bsm.vov, rho=sabr_bsm.rho, beta=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([25.631, 23.751, 21.924, 20.159, 18.465, 16.85 , 15.32 , 13.884,\n",
       "        12.543, 11.302, 10.159,  9.115,  8.165,  7.306,  6.532,  5.838,\n",
       "         5.217,  4.663,  4.17 ,  3.731,  3.342,  2.996,  2.689,  2.417,\n",
       "         2.175]),\n",
       " array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "        0., 0., 0., 0., 0., 0., 0., 0.]))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "price_hagan = sabr_bsm.price(strike, forward, texp)\n",
    "price_mc = sabr_bsm_mc.price(strike, forward, texp)\n",
    "\n",
    "# make sure the two prices are similar\n",
    "price_hagan, price_mc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat the same for beta = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'sigma': 20, 'intr': 0.0, 'divr': 0.0, 'is_fwd': False, 'vov': 0.5, 'rho': -0.25, 'beta': 0}\n"
     ]
    }
   ],
   "source": [
    "# instantiate mc model from the hagan model's parameters\n",
    "print(sabr_norm.__dict__)\n",
    "sabr_norm_mc = opt.sabr.ModelNormalMC(sabr_norm.sigma, vov=sabr_norm.vov, rho=sabr_norm.rho, beta=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([26.591, 24.753, 22.952, 21.192, 19.479, 17.818, 16.215, 14.675,\n",
       "        13.206, 11.811, 10.497,  9.268,  8.129,  7.08 ,  6.125,  5.262,\n",
       "         4.49 ,  3.807,  3.207,  2.686,  2.237,  1.854,  1.53 ,  1.257,\n",
       "         1.03 ]),\n",
       " array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "        0., 0., 0., 0., 0., 0., 0., 0.]))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "price_hagan = sabr_norm.price(strike, forward, texp)\n",
    "price_mc = sabr_norm_mc.price(strike, forward, texp)\n",
    "\n",
    "# make sure the two prices are similar\n",
    "price_hagan, price_mc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 3. Pricing under conditional MC method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "strike = np.linspace(75,125,num=25)\n",
    "forward = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'sigma': 0.2, 'intr': 0.0, 'divr': 0.0, 'is_fwd': False, 'vov': 0.5, 'rho': 0.25, 'beta': 1}\n"
     ]
    }
   ],
   "source": [
    "# instantiate mc model from the hagan model's parameters\n",
    "print(sabr_bsm.__dict__)\n",
    "sabr_bsm_cmc = opt.sabr.ModelBsmCondMC(sabr_bsm.sigma, vov=sabr_bsm.vov, rho=sabr_bsm.rho, beta=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([25.631, 23.751, 21.924, 20.159, 18.465, 16.85 , 15.32 , 13.884,\n",
       "        12.543, 11.302, 10.159,  9.115,  8.165,  7.306,  6.532,  5.838,\n",
       "         5.217,  4.663,  4.17 ,  3.731,  3.342,  2.996,  2.689,  2.417,\n",
       "         2.175]),\n",
       " array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "        0., 0., 0., 0., 0., 0., 0., 0.]))"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "price_hagan = sabr_bsm.price(strike, forward, texp)\n",
    "price_mc = sabr_bsm_cmc.price(strike, forward, texp)\n",
    "\n",
    "# make sure the two prices are similar\n",
    "price_hagan, price_mc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat the same for beta = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'sigma': 20, 'intr': 0.0, 'divr': 0.0, 'is_fwd': False, 'vov': 0.5, 'rho': -0.25, 'beta': 0}\n"
     ]
    }
   ],
   "source": [
    "# instantiate mc model from the hagan model's parameters\n",
    "print(sabr_norm.__dict__)\n",
    "sabr_norm_cmc = opt.sabr.ModelBsmCondMC(sabr_norm.sigma, vov=sabr_norm.vov, rho=sabr_norm.rho, beta=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([26.591, 24.753, 22.952, 21.192, 19.479, 17.818, 16.215, 14.675,\n",
       "        13.206, 11.811, 10.497,  9.268,  8.129,  7.08 ,  6.125,  5.262,\n",
       "         4.49 ,  3.807,  3.207,  2.686,  2.237,  1.854,  1.53 ,  1.257,\n",
       "         1.03 ]),\n",
       " array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "        0., 0., 0., 0., 0., 0., 0., 0.]))"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "price_hagan = sabr_norm.price(strike, forward, texp)\n",
    "price_mc = sabr_norm_cmc.price(strike, forward, texp)\n",
    "    \n",
    "# make sure the two prices are similar\n",
    "price_hagan, price_mc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Compare the MC variance between brute-force MC and conditional MC\n",
    "For this, you should not use the random number seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}