Equity Option IVOL.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Equity Option Implied Volatility Analytics with Python\n",
    "\n",
    "Modified notebook from Jason Strimpel's [talk](https://github.com/jasonstrimpel/PyData-Meetup).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Important jargon\n",
    "\n",
    "\n",
    "### Terms to know\n",
    "\n",
    "* __Call option__. Derivative contract which conveys the right (not the obligation) to the buyer of the option to purchase the underlying stock at a specific price on a certain date.\n",
    "\n",
    "* __Put option__. Derivative contract which conveys the right (not the obligation) to the buyer of the option to sell the underlying stock at a specific price on a certain date.\n",
    "\n",
    "* __In the money__. Condition where the stock price is greater than the strike price for a call option or lower than the strike price for a put option\n",
    "\n",
    "* __Out the money__. Opposite of in the money.\n",
    "\n",
    "* __At the money__. The stock price is equal to the strike price.\n",
    "\n",
    "* __Implied volatility__. Volatiltiy parameter that sets the pricing model equal to the market price of the option. (Much more on this later.)\n",
    "\n",
    "### Input parameters to price options\n",
    "\n",
    "Options priced with Black-Scholes have five input variables, four of which can be observed in the market and one of which is latent (unobservable, more on this later).\n",
    "\n",
    "* __Underlying stock price ($S$)__. Price of the underlying stock upon which the option derives its value.\n",
    "\n",
    "* __Strike price ($K$)__. The price at which the owner of the option can buy the underlying in the case of a call or sell the underlying in the case of a put.\n",
    "\n",
    "* __Time to expiration ($t$)__. Time (as a fraction of a 365 day year) until the expiration of the option.\n",
    "\n",
    "* __Risk free rate ($r$)__. So-called risk free rate (this is better described as the rate of interest you would earn in a riskless security over the holding period of the option)\n",
    "\n",
    "* __Volatility ($vol$)__. Volatility (or more usually implied volatility) is a latent value not actually observed in the market.\n",
    "\n",
    "Let's import our modules and declare these variables in Python for later use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numpy 1.14.0\n",
      "Matplotlib 2.1.2\n",
      "Scipy 1.0.0\n",
      "Pandas 0.22.0\n"
     ]
    }
   ],
   "source": [
    "# magic function to plot inline\n",
    "%matplotlib inline\n",
    "\n",
    "# python standard modules\n",
    "import time\n",
    "import datetime as dt\n",
    "from math import sqrt, pi\n",
    "\n",
    "\n",
    "# import numpy, pyplot and scipy\n",
    "import numpy as np\n",
    "import matplotlib as mat\n",
    "mat.style.use('ggplot')\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import scipy\n",
    "from scipy.stats import norm\n",
    "from scipy.optimize import brentq\n",
    "from scipy.interpolate import interp1d\n",
    "\n",
    "# import pandas\n",
    "import pandas\n",
    "# new in 0.17.1\n",
    "from pandas_datareader import data as pd\n",
    "\n",
    "\n",
    "print ('Numpy %s\\nMatplotlib %s\\nScipy %s\\nPandas %s' % \\\n",
    "    (np.__version__, mat.__version__, scipy.__version__, pandas.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jump to Yahoo fetch...\n",
    "- <a href=\"#fin\">Market Data Options Fetch</a><br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# underlying stock price\n",
    "S = 45.0\n",
    "\n",
    "# series of underlying stock prices to demonstrate a payoff profile\n",
    "S_ = np.arange(40.0, 50.0, 0.01)\n",
    "\n",
    "# strike price\n",
    "K = 45.0\n",
    "\n",
    "# time to expiration (you'll see this as T-t in the equation)\n",
    "t = 164.0/365.0\n",
    "\n",
    "# risk free rate (there's nuance to this which we'll describe later)\n",
    "r = 0.02\n",
    "\n",
    "# volatility (latent variable which is the topic of this talk)\n",
    "vol = 0.25"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What are options?\n",
    "\n",
    "Options are standardized derivatives contracts that convey the buyer the right (but not obligation) to buy the underlying security in the case of a call option or sell the underlying security in the case of a put option at a given price before a given date (for American style options) or on a given date (European style options).\n",
    "\n",
    "We will focus on plain vanilla, european style, equity options so we can use the famous Black-Scholes pricing formula. More on Black-Scholes a bit later.\n",
    "\n",
    "Most derivatives have a payoff function which describes the value at the end of the life of the contract (expiration). For a call option, the payoff can be expressed simply as:\n",
    "\n",
    "> $$C(S, K) = max(S-K, 0)$$\n",
    "\n",
    "While the payoff for a put option can be expressed as:\n",
    "\n",
    "> $$P(K, S) = max(K-S, 0)$$\n",
    "\n",
    "Where $S$ and $K$ are defined above. In this case, $K$ remains fixed during the life of a contract while the underlying stock price, $S$, fluctuates..\n",
    "\n",
    "Across a range of potential $S$ values, we can form what is commonly known as the payoff (or risk profile or PnL chart) for an option at expiration.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7a6fe200f0b8>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEaCAYAAAAxEb1YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XecXFXdx/HPWTZBeltagCSUAIIUQamCCFJEHkPRoySAASSAdB76g0iVSEelhYChBTiIVEFFJCCIlCBICU0SWoCwlPSQZOc8f5wzm9lhdudumblzZ77v12teu3Pnlt+5d+Y3Z8698xvjvUdERLKjKe0ARESke5S4RUQyRolbRCRjlLhFRDJGiVtEJGOUuEVEMkaJO2XGmMnGmNML7o83xoxJMZ6xxpi/pbX9tBhj1jDGPGyMmWWM8XHa0saYu4wx040x3hgzON0o+06ln2fGmBHGmAWVWn+jU+LuIWPMCsaYC4wxrxlj5hpjphpjHjPGHGCMaU47vnKMMfvlE1SRY4AfVTueGnAasBKwKbBqnHY4sDWwbZz2bk9WbIw53RgzuQ9i7Et7A8f3diXGmNXjm9oORQ/dDqzW2/VLaTWfYGqRMWZ14AlgAXAG8G9gPrANcALwH+D51ALsBe/9tLRjSMkQ4Gnv/RtF01723r+YUkwV473/tKvHjTH9vffzerH+OcCcni4vZXjvdevmDbgP+BBYpsRj/YAl4v87A+OBT4FpwKPAFkXzTwZOL7g/HhhTZvvrAX8CZsbbfcA6BY+PILypfBd4GZgLPA1sFh/fAfBFt7HxsbHA3wrWZQhvRm8B84D/AseWaMPZwOWxrR8BFwGLdNGGwXG7+wMPE17kk4DhRfOdB0wEZhN6vFfn9zuwFDADGFZi3Tlgh4L5rgE+jvviWWCXgvm/tC9imwqnje+kHQa4Nu6XOXE//QpYtOBYFK//zC72yzrAncDnwGfAX4GNCh6/Ksa2bMG03wNvAEsVHI/zgDHAdKAV+DXQ1NnzLN6/DjgH+AD4OE4fBjxFeP62Ep5363ax7yYXPgeL2rY7MAH4ApgKXEl8rRQ+94CRwNsx9nuAFdN+zdfaLfUAsnYDlgfaKEi2Xcy7F2HYYV1gw/hC+hRYoWCeyXQjcQOLxSf1w8Dm8fYI8CbQP84zgpC4ngO+DWwM3B9fkIsD/YEj4gttlXjLJ8OxdEzcRxAS0khCD/QwQvI7uKgNnwGnxHl+THjjOLCLdgyO258CDCe8GZ0b4/5GwXynA9vF+XcCXgVuKHj8GuCRonWfA7wOmHj/jhjjrsBXCW8w84D14+OrAP8EbsnvC2BFwsf9x+K05TtpR1OMe8sY4w/ifj6r4HiNIrzp5Pf1kp2sa2VCh+AqYKO4T34LfEJMXsBXCJ/o7oj3hxES4eZFx2M64c10PcKb4yzg+M6eZ/H+DMIb4wbENwvgQGAPYG3g68C9hDeJ/HPt6/E47h3blo9zBAWJm/AcXABcGo/B94B3gJsK5hlLeIO4Ffga4RPs24XHW7e4r9IOIGs3YIv8E7UHyzYREtzwgmmT6V7iPpjQ+2wpmLYyIbkeEO+PiDHuVDDPcoTe+c/i/f0AX2L9Y+mYuN8FLiia51LgraI23Fs0z5+BW7tox+AY4zlF0/8J3NzFcnvFRNUU728W1zMk3l8kxnxivL9OfHz3ovU8B1zf1X4v3hfdOM7HAW8U3D+d2BMts9yZwL+KphmKPuXExDcLOJ+QoI8rWmYy8I+iab8C3uusvfH+6xT0yjuJcfm4P7eN91eP93comm8EHRP3TYShqMJ5hhLeqAcV7O+PiZ9W4rRTgA+6ewzq/aaTk91n4l9fdkZj1jTG3GSMedMYM53wIlsGGNSL7W8IvOK9b81P8N5/BLwWHyv0ZME8nxGGHDZIuiFjzNKEF+ZjRQ89Cgw2xixeMK14TP99whtKOU8W3X+iMEZjzN7xpO8UY8xMQq+4P6F3h/f+OcLQx8/iIt+L270h3s+vq7gNj/Hl/dUjxphDjDFPGWM+ijGeT8+O8TeBzY0xM/M3Qi94MOGTDADe+4mE4atTgMeBy0qsq9R+XS0e085M8N7nCicYYzaNV9ZMMsbMIPSSofvt25DSzyNDx+fkRO/9FwX3kz6PGooSd/e9QeglJHnR3w8MJAw3bEW4YmEqIfH0Rqk3DdPJ9OJ5+mJ7pdZTfCLL07PnV/u6jTFbEoY5HiP0tDcjDNVAx314NTDCGNOPkMDv9t5PTbCdsm++ZYM15kfAFYRhld0JQwdnE851dFcTYQhs06LbeoTeeKHtCUN2AwnDJ2VDTTDPrA4LhDfmvxL200GET5vfjPd78hzubH8XTi/1POrp87ZuKXF3kw9n4x8EjjTGLFP8uDGmnzFmCWPMCoSexCjv/V+8968QxoZX6mUILwMbGmNaCra5MmEc/eWiebcqmGdZYH1CrxviC8QYs0hnG/LeTwfeI4yTF9oemOS9n93DNpSMMdq6IMZvAa3e+9O99095718nfAIodhsheR0KfJ9wsjAvv0+2L1pmO768v3pie+Df3vtLvPcTfLgqZXDRPPMIQzjlPEvoELzvvX+z6PZxfiZjzMHAnoTjsjhh6KpYqf06JR7TpL5KGOv/P+/9I7GnvxwdE2k+0ZZr38t8+Xn0bUJifqUbMQlK3D31c8LlfxOMMcOMMRsYY9YxxuxHePENIYxlfwwcYoxZ1xizNeGkS28vkRoX13u7MWYzY8zmhMT1PqHXl+eBC4wx2xtjNgJuJPSoxsXHJ8W/PzDGrGiMWbKT7Z0PHBWHA4YYYw4lXN/8q162I+/guA/XNcacTUgw+Y/+rwErGmMONsasZYw5gLDvO/DezwJuBi4mfJT/W8Fj/yX02q80xuxqjFnfGHM54eTXhX0Q/2vARsaYocaYtY0xxxBO1BWaBKxijNnaGNNSNMRU6HeEBHi3MWY7Y8xgY8y3jDHnGWO2ATDGrEc4uXqc9/4JYF/CPtynaF2bGmPOjPt1GOH6/FIJvitvE84nHBXbtlPcdmEPuZVw7mQXY8wqxpjlOlnXhcBmxphL4jHYjXDi9Rbv/TudLCOdSXuQPas3Qk/kYsIJnbmEIZBHCSf9muM83wZeiI+/BuxDuPrjzIL1TKZnlwM+wMLLAe+n9OWAuxB6r18Az1BwtUac7zLCpXs5ur4c8ERC8plPuNyt1OWApxdNG0Mnl9DFxwez8HLA8XEfTQb2L5rvnBjjrNjmfeNyg4vm2yROP7XEtpZm4eWAX1B0OWBn+714X3TSjn5x3Z8SzmGMA46k4MRvnGdcnMfT9eWAgwjj+PlY3ya8Ka0JLEr4zsCdRcucSugoDCo4HucRLhOcHrd7IQWXZxa3t7PnHfBDwvDg3Ljtb8fn1oiCeQ4oeH5MLnwOFq2r8HLAjwlXz3zpcsCiZUqeRG/0W/5yKakjxpgRhBdhzX7BKn59fBKwnff+8T5Y3+7A3cBA7/2HvV1flsVvaY7x3p+bdixSGTX7whZJIg47DCR8g3VcoydtaQwa45asOwl4iTDcc1LKsYhUhYZKREQyRj1uEZGMqdQYt7rxIiLdl+jLRhU7OTllypQeLdfS0kJra2v5GeuI2twY1Ob615v2DhgwIPG8GioREckYJW4RkYxR4hYRyRglbhGRjFHiFhHJmERXlVhrJxMKurcBC5xz36hkUCIi0rnuXA74Hedc41zXIyJSozRUIiLSB/zrLzPr3tuoRhmRRLVKrLWTCPV+PXCNc250iXlGEn4JHOfc5vPmFf8CUTLNzc0sWLCgR8tmldrcGNTm+tX22Sd8+r8jMIstzgoXj8V8ZbFur6N///6Q8JuTSRP3AOfcFGvtSsBDwFHOueIf/izk9c3J5NTmxqA21yff1kbu0jNg0mssf8F1TFviS79omEj85mSixJ1oqMQ5NyX+nQrcRfjRUBGRhufvuRleexGz3xH0G7R2VbZZNnFba5ew1i6V/5/wc1gvVTowEZFa559/Cv/gnZjtd6Vp6+9UbbtJripZGbjLWpuff5xz7s8VjUpEpMb5jz8kd/1lMHBtzE8Oqeq2yyZu59xbhB9iFRERwM+fR+7qUWAMTYedjOnXv6rb1+WAIiLd5G8dDe+8RdPBx2FWXKXq21fiFhHphtwTD+P/8VfM7j/CbPzNVGJQ4hYRSci/Owl/y1Ww/saYocNSi0OJW0QkAT97VhjXXmJJmg45AdO0SGqxKHGLiJThvSc39nJo/YimQ0/CLL1sqvEocYuIlOEfuhv+/S/MPiMw62yQdjhK3CIiXfGvv4y/8wbYbBvMzkPTDgdQ4hYR6ZSf9hm50RdCyyo0jTgaYxKVEqk4JW4RkRJ8Wxu5ay+COTNpOvwUzGKLpx1SOyVuEZES2otHDf85ZvXBaYfTgRK3iEiRDsWjttkx7XC+RIlbRKRAmsWjklLiFhGJ0i4elZQSt4hI1F486qB0ikclpcQtIgLk/hmLR33vh5hN0ikelZQSt4g0PP/eJPzNV8F6G2GGDk87nLKUuEWkofnZs8hdNQoWX5KmkSdgFkmveFRSStwi0rC89+Ru+E1B8ajl0g4pESVuEWlY/qG74bknQ/GoIekXj0pKiVtEGlItFo9KSolbRBpOrRaPSkqJW0QaSi0Xj0pKiVtEGoq/55ZYPOrwmiselZQSt4g0DP/C0/gH/4DZbheattkp7XB6TIlbRBpCKB51aSgete/ItMPpFSVuEal7oXjUrwFqunhUUkrcIlL3QvGo/9J00PE1XTwqKSVuEalrWSoelZQSt4jUrawVj0pKiVtE6lIWi0cl1Zx0RmvtIsCzwPvOuT0qF5KISO90KB51wq8yUzwqqe70uI8BJlYqEBGRvuIfuicWj/pppopHJZUocVtrVwe+D4ypbDgiIr3j33gFf+dY2GxrzM57ph1ORSQdKrkMOAlYqrMZrLUjgZEAzjlaWlp6FlBzc4+XzSq1uTGozZXX9vmnfHrtRSyy8gCWP/4smpZYsmrbhuq1t2zittbuAUx1zk2w1u7Q2XzOudHA6HjXt7a29iiglpYWerpsVqnNjUFtrizf1kbu0jNg1nSajv4Fn86ZC3PmVmXbeb1p74ABAxLPm2SoZFvgB9baycBtwI7W2pt7FJmISIV0LB61ZtrhVFTZHrdz7lTgVIDY4z7BObdfheMSEUmsXopHJaXruEUk0+qpeFRSia/jBnDOjQfGVyQSEZFuqrfiUUmpxy0imbWweNRxdVE8KiklbhHJpIXFo/bBbLJF2uFUlRK3iGSOf28y/pZ88ajGu1ZCiVtEMqW9eNRi9Vc8KiklbhHJjIXFoz6kaeSJdVc8KiklbhHJjA7Fo9bdMO1wUqPELSKZ0AjFo5JS4haRmuenf0bumgugZWWafno0xpi0Q0qVEreI1DTf1kZu9EUweyZNh5+CWXyJtENKnRK3iNQ0f++4hikelZQSt4jULP/CM/gH7gjFo7at/+JRSSlxi0hNCsWjLoGBazVM8aiklLhFpOZ0LB51SsMUj0pKiVtEao6/7dqGLB6VlBK3iNSU3D//jn/sLw1ZPCopJW4RqRmheNSVDVs8KiklbhGpCSoelZwSt4ikTsWjukeJW0RS5/92bygetXdjF49KSolbRFLl34zFo76+FWaXxi4elZQSt4ikxk//PBSPWmElmkYc0/DFo5JS4haRVPhcG7lrL4JZM8OXbFQ8KjElbhFJhb9nHLz6n1A8ag0Vj+oOJW4RqToVj+odJW4RqSoVj+o9JW4RqZr24lFexaN6Q4lbRKpmYfGoY1U8qheUuEWkKtqLR+22D2bTLdMOJ9OUuEWk4joUj9pTxaN6S4lbRCrKz5m9sHjUISoe1ReUuEWkYrz35MYWFI9aRsWj+kJzuRmstV8BHgMWjfP/wTn3y0oHJiLZN/u+2+G5f2J+eKCKR/WhJD3uL4AdnXObAJsCu1lrt6psWCKSdf7NV5h54xUqHlUBZXvczjkPzIx3+8Wbr2RQIpJt+eJRi6y0Kl7Fo/qc8b58DrbWLgJMANYBrnDOnVxinpHASADn3Obz5s3rUUDNzc0sWLCgR8tmldrcGBqlzb6tjc/POpZ5r73IShddj1ljrbRDqpreHOP+/fsDJHqHS5S486y1ywJ3AUc5517qYlY/ZcqUxOst1NLSQmtra4+WzSq1uTE0Sptzd90U6pCMOJqVhv6kIdqc15tjPGDAAEiYuLt1VYlz7nNgPLBbt6MSkbrn/xOLR31rZ5q2/W7a4dStsonbWrti7GljrV0M+C7waqUDE5Fs8a0fkbvuUlhjTRWPqrCyJyeBVYEb4jh3E+Ccc/dXNiwRyZKFxaM8TYefium/aNoh1bUkV5X8B/h6FWIRkYzyt42Bt9+k6Yj/U/GoKtA3J0WkV3JPPoJ/7M8qHlVFStwi0mP+vcn4m69Q8agqU+IWkR5ZWDxqCRWPqjIlbhHpNhWPSpcSt4h0m//bvaF41N4HYNb9WtrhNBwlbhHpFv/mK/g7x8KmW2F22SvtcBqSEreIJJYvHsXyK9J04NEqHpUSJW4RScTn2siNuRhmzQxfsll8ybRDalhK3CKSiL/nVpj4Amb4YZg11kw7nIamxC0iZYXiUU7Fo2qEEreIdEnFo2qPEreIdKpD8ajDTlHxqBqhxC0inWovHnXQMZiVVk07HImUuEWkpPbiUbvujdlUvw9eS5S4ReRL2otHrfs1zF77px2OFFHiFpEO/JzZYVx7sSVCHRIVj6o5Stwi0s57T+6G38DHH6h4VA1T4haRdv7he2GCikfVOiVuEQFi8ag/jFXxqAxQ4hYRFY/KGCVukQbXoXjUYaeoeFQGKHGLNLj24lHDDsUMXCvtcCQBJW6RBtZePGrb79L0rZ3TDkcSUuIWaVAdikcNOzTtcKQblLhFGpCfP1/FozJMiVukAfnbr1XxqAxT4hZpMLl/PYJ/VMWjskyJW6SB+Pffxt+k4lFZp8Qt0iD8nNnkrhql4lF1QIlbpAF0KB51iIpHZZ0St0gDaC8etdf+mPVUPCrrmsvNYK1dA7gRWAXIAaOdc5dXOjAR6Rv+zYmxeNSWmF33Tjsc6QNJetwLgP91zn0V2Ao4wlq7QWXDEpG+0LF41DEqHlUnyiZu59wHzrnn4v8zgInAapUOTER6Z2HxqBkqHlVnyg6VFLLWDga+DjxV4rGRwEgA5xwtLS09C6i5ucfLZpXa3Biq3eaZ40Yza+ILLH3EaSy22RZV226hRjvO1Wqv8d4nmtFauyTwKHCec+6PZWb3U6ZM6VFALS0ttLa29mjZrFKbG0M12+xffJbcb84OxaNGHF2VbZbSaMe5N+0dMGAAQKKxrERXlVhr+wF3ArckSNoikiLf+hG5MZeoeFQdK5u4rbUGuA6Y6Jy7pPIhiUhPqXhUY0gyxr0tsD/worX2+TjtNOfcA5ULS0R6or141M9PU/GoOlY2cTvnHifhuIuIpGdh8ai9MF9X8ah6pm9OitSBUDzqSlh3Q8xeB6QdjlSYErdIxvm5s8ldPQq+slioQ6LiUXVPiVskw7z3+LG/hakf0DTyJMyyy6cdklSBErdIhvmH78NPeELFoxqMErdIRoXiUb9X8agGpMQtkkEqHtXYlLhFMqa9eNTM6TQddrKKRzUgJW6RjPH33goTX8AMOxQzcO20w5EUKHGLZIh/8Vn8nxxm251o2m6XtMORlChxi2SE/2RqKB61+pqYYYelHY6kSIlbJAP8/PnhF9p9jqbDT1bxqAanxC2SAd6NCcWjDjwWs9KAtMORlClxi9S43L/G48c/qOJR0k6JW6SG+fffwd90hYpHSQdK3CI1KhSPOl/Fo+RLlLhFalB78aiPPqBp5IkqHiUdKHGL1KCOxaM2SjscqTFK3CI1pr141CZbYHZT8Sj5MiVukRriZ0xbWDzqoGNVPEpKUuIWqRE+10bu2otUPErKUuIWqRH+vttUPEoSUeIWqQH+xQn4+29X8ShJRIlbJGX+k6nkrrsEVh+s4lGSiBK3SIrai0fl2mg6/BQVj5JElLhFUtRePGrEMSoeJYkpcYukpL141C57YTbbOu1wJEOUuEVS0F48asgGmL1VPEq6R4lbpMo6FI8aeZKKR0m3KXGLVJH3Hn/D71Q8SnpFiVukivzf78c/+7iKR0mvNJebwVp7PbAHMNU597XKhyRSn+a9+iL+jutVPEp6LUmPeyywW4XjEKlrfsY0pl30C1iuRcWjpNfKJm7n3GPAp1WIRaQu5YtH5aZ/Hr5ko+JR0ktlh0qSstaOBEYCOOdoaWnpWUDNzT1eNqvU5vo2c9y1zJr4AsseeRqLbrZl2uFUVSMdZ6hee/sscTvnRgOj413f2trao/W0tLTQ02WzSm2uX/7FCeTu+D1mm53ov+P3G6LNhRrlOOf1pr0DBiT/5qyuKhGpkOLiURrXlr6ixC1SAX7+fHJX/3ph8ahFVTxK+k7ZxG2tvRV4EljPWvuetfbgyoclkm3eXQeT31DxKKmIsmPczrl9qxGISL3IPfUofvwDKh4lFaOhEpE+5N9/B3/j71Q8SipKiVukj6h4lFSLErdIH1DxKKkmJW6RPrCweNR+Kh4lFafELdJL/r+vLiwetauKR0nlKXGL9IKfMY3cNReE4lEHHotp0ktKKk/PMpEe8rk2cmMuhhnTwpdsllDxKKkOJW6RHvL33Q6vPI8Zdihm4NpphyMNRIlbpAf8SxPwf7ods81OmG/tnHY40mCUuEW6yX8yldyYS2C1QSoeJalQ4hbphg7Fow5T8ShJhxK3SDcsLB51NGZlFY+SdChxiyS0sHjUnpjNtkk7HGlgStwiCbQXj1pnA8xeKh4l6VLiFimjQ/GoQ0/ENPfZL/6J9IgSt0gXvPf4G68oKB61QtohiShxi3TF//1P+Gf+oeJRUlOUuEU6oeJRUquUuEVKWFg8agUVj5Kao2ejSJEOxaMOU/EoqT1K3CJF2otH7TsSM0jFo6T2KHGLFGgvHrX1jpjtdkk7HJGSlLhFog7Fo4YfruJRUrOUuEWIxaOuuUDFoyQTlLhFAH/HdTDpdRWPkkxQ4paGl3vqUfwjKh4l2aHELQ3NT1HxKMkeJW5pWH7ubHJXjYJFv6LiUZIpStzSkBYWj5qi4lGSOUrc0pDai0ftORyz/sZphyPSLUrc0nDai0dt/E3MbvukHY5ItyUa1LPW7gZcDiwCjHHOjapoVCIV0qF41EHHqXiUZFLZZ621dhHgCuB7wAbAvtbaDSodmEhfC8WjLlHxKMm8JD3uLYA3nXNvAVhrbwOGAq/0dTBt5x5Pa66NtgUL+nrVNa21uVltroYF8+HjDzH7H6HiUZJpSRL3asC7BfffA7YsnslaOxIYCeCco6WlpdvBTBu0FrQtoDnnu71slpkmozZXSb/v/5DFf7BvKnVImpube/S6yLJGa3O12pskcZd6hn/pFeecGw2Mzj/e2tra/Wj2P5KWlhZ6tGyGqc3VswCY88knVd8u6Dg3gt60d8CA5KUWkpyZeQ9Yo+D+6sCUbsYkIiJ9JEmP+xlgiLV2TeB94CfAsIpGJSIinSrb43bOLQCOBP4CTAyT3MuVDkxEREpLdB23c+4B4IEKxyIiIgno2wciIhmjxC0ikjFK3CIiGaPELSKSMcb7inx7rbG+Bigi0jcSfaW3Uj1u09ObtXZCb5bP4k1tboyb2lz/tz5obyIaKhERyRglbhGRjKnFxD26/Cx1R21uDGpz/atKeyt1clJERCqkFnvcIiLSBSVuEZGMSVRkqpLib1o+C7zvnNsjlo+9DVgeeA7Y3zk3L80Y+1qJNt8CfAOYDzwNHOqcm59mjH2tuM0F038LHOicq7sfgCxxnA1wLvAjoA24yjn3mzRj7Gsl2rwTcCGhkzgTGOGcezPNGPuStXYyMINwPBc4575hrV0euB0YDEwGrHPus77cbi30uI8hlIvN+zVwqXNuCPAZcHAqUVVWcZtvAdYHNgIWA36WRlAVVtxmrLXfAJZNJ5yqKG7zCMKPkqzvnPsqoYNSb4rbfBUw3Dm3KTAOOD2VqCrrO865TZ1z34j3TwEejjns4Xi/T6WauK21qwPfB8bE+wbYEfhDnOUGYM90oquM4jZDKJvrnPPOOU/oca+eVnyVUKrNsWd2IXBSWnFVUqk2A4cDZzvncgDOualpxFYpnbTZA0vH/5ehMX49ayghd0GFcljaPe7LCC/cXLy/AvB5/PEGCD+btloagVVQcZvbWWv7AfsDf652UBVWqs1HAvc65z5IJ6SKK9XmtYEfW2uftdY+aK0dkk5oFVOqzT8DHrDWvkd4bo9KI7AK8sBfrbUT4g+mA6ycf17Hvyv19UZTS9zW2j2Aqc65CQWTS33ls26uV+ykzYWuBB5zzv2jimFVVKk2W2sHEMZ5f5taYBXUxXFeFJgbP1JfC1xf9eAqpIs2Hwfs7pxbHfg9cEnVg6usbZ1zmwHfA46w1m5fjY2m2ePeFvhBHNy/jTBEchmwrLU2f9K03n6Y+EttttbeDGCt/SWwInB8euFVRKnj/DKwDvBmnL64tbZuTljR+XF+D7gzznMXsHE64VVEqTb/CdjEOfdUnOd2YJuU4qsI59yU+Hcq4ZhuAXxkrV0VIP7t8yGx1K4qcc6dCpwKYK3dATjBOTfcWnsH8EPCwf8pcE9aMfa1Ttq8n7X2Z8CuwE758c960Umb9yicx1o70zm3TgrhVUQXx3kU4Y3reuDbwOupBdnHSrWZMLb7obV2Xefc68DOFJ2gzjJr7RJAk3NuRvx/F+Bs4F5C7hpFhXJY2mPcpZwMHB97YCsA16UcTzVcDawMPGmtfd5ae0baAUlFjAL2sda+CJxPfV491C6eqzoEuNNa+wJhjPvEdKPqUysDj8e2PQ38yTn3Z8Jx3tla+wbhzarPx/X1lXcRkYypxR63iIh0QYlbRCRjlLhFRDJGiVtEJGOUuEVEMkaJW0qy1o6w1j7ei+XHWmvPTTDfQGvtzFi7JHOstePjdfiVWPd21trXKrFuybbUy7pK37DWemBIYclMa+2ZwDrOuf1SC6wM59w7QFVLusYviNwcv4Zds2Lpg/UquQ1rbX/CNeU/JlRq/Bi42zl3XCW3K72jHrf0uaz2nmtJQdmHSjuVUAt+C2DeTMGFAAAGR0lEQVQp4DvAv6u0bekh9bgbRL6XCVxK+HZqG3Cac+738fEVCEWAdgBeBf5StPz6hKJQmxN6Zb9wzrn42FhgDjCI8FXuoUXLvgSc6py7L97vB3wAfBf4HJgE9HPOLbDWjgf+Qfhq+MbAk8Aw51xrXPYA4BxCL/0yQr32nznn/laizbsDFxFqYE+Pbb8KeBBY1Fo7M866LvAJoRa8jdMccLJz7ou4rqHAWcBasf1HxG/JFW5v1bjfbnTOXVQinsnANYRvEK4K3A0c7pybW3B8fksozPSQtfY6Cj4ZWGvXAC4HtiN0um51zh0ZHzuI8K3EVQjf4hvpnHu7OIYSvgncla+5QSj8PznBcpIi9bgbyyqEmsirERLeFdba5eJjVwBzCQnloHgD2msyPEQohL8SsC9wpbV2w4J1DwPOI/TaisfGbwQKh2t2Bz5wzj3fSZzDgAPjtvoT6l5grd2AUEFxeIwz35bOXEf4NaGlgK8Bf3fOzSJUcpvinFsy3qYA/wdsBWwKbELogZ4et7tFbMOJhOGE7SlKbtbawcCjwO9KJe0Cwwl1adYmvGEU/rDAKoRffhoEjCxcKH6KuR94m/DLKqsRf4jBWrsncBqwN6FQ2T+AW7uIodC/CCUmfm6t3SjWxJcapx53Y5lPKOS/gFAjeSawnrX2GWAfYKOY2F6y1t5ASFAAewCT871z4Dlr7Z2EYmAvx2n3OOeeiP/PtTbfcQVCT/IX1tqlnXPTCT3Om7qI8/exKBHWWgf8IE7/IXCfc+7x+NgZwNFl2ruBtfaF+NNRXf181HDgqPyPG1hrzyL0jn9BeJO73jn3UJz3/aJlNyAk4FOdc+US5u+cc+/GbZxH6GHnk3cO+GVBL79wuS2AAcCJBfXq82+QhwLnO+cmxuV+BZxmrR2UoNd9PmG/DCd8IvnEWnuqc+6GrheTNKnHXT/agH5F0/oRklfeJwUveoDZhCGHFQlv4u8WPFb4gh8EbGmt/Tx/I7zQVymYp3DZDmKP9glCgaVlCT3eW7poy4clYoSQuNq345ybTRji6Mw+hN7929baR621W3cx7wA6tvntOA3CUMt/u1h2OCGZ/6GLefKK9/GAgvsfO+fmdrLcGsDbRccvbxBwecGx+ZRQ277sj5A459qcc1c457YlfJo4D7jeWvvVBG2RlKjHXT/eIXyELiybuSbJSod+DCwgJIdX47SBBY+/CzzqnNu5i3WUq1Z2A6EaXjPwpHOuuNeaxAcUXGVhrV2MUEGyJOfcM8DQOKZ+JGHceo1OYp1CSID5TxADWVgL/l3C0EZnzgR2A8ZZa3/inGvrYt41Cv4v3AadxJX3LjDQWttcInm/C5znnOvqzbAs59wcwvDZWYRPEXVTgrXeKHHXj9uB02PJ0CmEk3v/A3TVywRCr8ta+0fgzHiSazChjvDkOMv9wChr7f4s/IHbTYGZ+Y/nCdxNGJ9eGbgg4TLF/gD8y1q7DeGXxM+i9K8m5S9z+xFwv3NumrV2OuFTCcBHwArW2mWcc9PitFsJ++8ZQgI9gzDEA2Gs/K/W2vuBRwjj60s55/JvcvPjtu4GbrLW7tdFXfUj4npmE8alb0/Y9qcJb1yj4o9utAGbx+Gpq4FzrLXPO+dettYuA+zinLsj7ovxwHjn3Jkl9tOxwPPAU7EdwwnnKXRlSQ3TUEn9OBv4J2Hc8zNCchzunHsp4fJHEoYkPgTGEq4wAcA5N4NQJP4nhDeFDwlXYCyaNLjYm7uT8Cngj0mXK1rHy8BRhDePD4AZhF8X+aKTRfYHJsekfRjxBGlMuLcCb8XhhQHAuYQ3g/8ALwLPxWk4554mnCy9FJhGOAk5qCi2eYSTgysRhho6e22NA/4KvBVvZb+kFNffRngjXofw6eo9wrXXOOfuIhyP22JbXyIMR+WtQRiqKmUOcDHhmLYCRwD7OOfeShKXpEP1uKVq4snEdfvqC0HW2iUJlxMOcc5N6ot1VlK8HLDkpYsV3ObqwB3OubKfvCQ7NFQiVWGtXZ5wdcb+vVzP/wAPE4ZILiL0jif3Nr565Zx7jwTDZZItGiqRirPWHkI4gfagc+6xXq5uKGG4ZgowBPiJc04fG6WhaKhERCRj1OMWEckYJW4RkYxR4hYRyRglbhGRjFHiFhHJmP8H6uFAKBPnYoUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7a6fe98e7e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEaCAYAAAAxEb1YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYU9X9x/H3NwzgBooMiiOiVKBKtWpd6i7WfSst2qOCC26IW9W60vbXalsrWq3VakUUBXePu0WtW8WldbfuWouKG7ggKiDrkPP749yxaczMZJbMnSSf1/PkgST35n5Okvnm5Nybcy2EgIiIlI9M2gFERKRlVLhFRMqMCreISJlR4RYRKTMq3CIiZUaFW0SkzKhwCwBmNtTMgpn1SztLRzOz35vZx0n7RyW3HWdmH5hZ1szOSDdh+zGzUWZWX+JtBDM7oJTbqHYq3G1gZpOSN2kws3oze9fMxptZ7xY+Tn1DwegIjWzvn8BqwIyOytEZmNn3gbHAaGL7bzKzOuBPwNnA6sB5rXzsfsl7Y2g7xW0PNxHb1GZm9qCZTSpw12rALe2xDSmsJu0AFeAxwBGfy42BK4A1gD3SDNVSIYTFwEdp50jBICAbQriz4QYz25jYqbkrhDAztWQlEEJYACxo7H4zM6AmhLCkDduoxvdRxwoh6NLKCzAJeDDvtl8AS4FlgbWAAGydt8w04Izk/9OTZb6+NLG9rsA44ENgMfAaMCJvmQAcD9wKfEXsQf8s5/6C2wOGJtf75Sy7OfAo8Q/9c+B6YJWc+89I2jIMeCPZ3sPA2s08b1OBK5O2zALmED/wls1ZZqdkudnAl8AjwGY5908G7i/w2A8Dk3KuH5w8T4uAD4DfEQtTw+sX8i5nFLhtrUbaMQJ4Ksk3C7gbGJz3WuRepjfxnNQk234HWAi8ChyZc/+WwBLgxzm3bQ/UA7vnvR4jgLeTx3kQGJCzziigPv968lj/Ir6v9gQGALcl75/5wMvAgXnv/fz2Dc1p9wE5y64G3Ah8kbyXpgKb5Nw/NFlnJ+L7bX7ymu2S9t94Z72kHqCcLxQu3D9L3oQ9KK5w90n+cI4H+gJ9m9jeH4DPgJ8Ag4GfA1lgh5xlArHYHZcsc3zy+MOb2h55hTu5bw6xWK8PbA28BDyWs60ziMX6b8RvGxskf/yPNPO8TU0e+3JgXWAv4BPgopxlfpzTzu8QC/tsoHdy/xZJ23OL0trJbVsl1/cgfoiOTR5nX+IH0G+T+1fMeX76JpcVgOHJc7FRcluXRtpxCLHIrZ0sexfwH6Bbcv9GyeMMTx6nTzPvpZeAnYlFc19ioTssZ5lfJM9Bf6CW+AF+foHX43Fg0+TyFPAiYMkyo/hm4c4CzwA/AL6VvEfWB44Bvpu077jkedo+57l7lDj00vDcNbT768INWJLhBeJ7aP1knc+B2rz33ovArsRvQVcn7V8p7b/zznhJPUA5X8gr3MAQ4C3gyeT6WjRTuJPr9cCoZra1HLHXeHTe7bcDf8+5HoBr8pa5Hni8qe3xzcL9W2IPtVvOMhsky2ybXD8jeaw+OcvslxSCZZpoy1Riz79Lzm2jk/Yt38g6meSPfWTObS8Bv8u5fjbwas71xwCf9zjHE3t9DUVmFDmFrNBz0YL3w8rJeg0fHP3I6Yk2sd6A5DlbJ+/2XwEv5D0HDybtuodYbLvm3H9Gsr2BObcNTm7bsVB7k+sB2KaI9t0JXJ5z/UFyvt3kvQcbCvcOyfUhOfd3B2YCv8p7vofnLNM3uU297gIX7Zxsu6FmNs/MFgCvEL+ijijBdgYC3Yi9nFyPEHukuZ7Iu/4P4odKS3yH+AG0uOGGEMKLxGGB3O3NCCF8mnP9Q2Iva5VmHv/pEMLSvIzdiL07zGyAmV1jZtPMbA6xh74isGbOOpcBh5hZFzOrIRahy/PaUOj5WqZhO21hZhua2e1m9o6ZzQXeS+5as6n1CtiE+Jw9m7yX5pnZPOI3qkENC4UQssCBxF7rtsB+4Ztj0Z+GEKblrPMmcRinudf/mby2LWdm48zsVTObneTZvRVt+w7wWQjhtZxMi4i98Pz37Qs5y3xE/La0agu3VxW0c7LtniKOo9YDM5M3ZYNs8q/lrdO1DdsLedetwG358rff2m0Vun1xI/e1tFOQn3EKseAcA7yfbOdxYnFvcA1wDnFIJAP0In7Fbixr7naae86aDmu2HHB/kulQ/rtj99W8jMVoeK62JI7v5srPuSFxOCcQd4K/VUzcZu5fGkJYmHfbH4j7Lk7iv/svzid+eLZUoee60Ps2/70EOvKtID0pbbcghDAthDA9r2gDNPRE6xpuMLNV+ObhWIuBLs1sZxpxKGG7vNu3JRaLXJvnXd8CeL2F23sV2MLMvi5CZrYB8Q83f3utsamZ5WbYIsn1VnI45RBgXAjhvqS3tpC8XnwIYQ5xp9cRyeXWEMLsvDYUer4WEL8ZtcW6xLHgX4QQHg4hvE784Mgtkg2FqLnn+rnk3/7Jeyn38nVhNrO+xJ2yZwMXAdcWOPS0j5mtnbPOYKA3//v6F2Nb4LoQwk3JN623icMuuYp9H9Wa2dc9fjPrDmxG+7yPqpIKdwmFeOjVP4BTzWyD5DCzq4kFONc7wPZmVmdmtY081nziH+tvzewnZjbIzH5O7BX9Pm/xPc3s2GSZ44g7ui5oyfaAi4GewCQzW8/Mtib2cB8PITxW7HPQhN7AJWa2rpntQRxTvzyE8BVxLPtT4AgzG2xmWwA3UPgwtsuA3YBdgAl5950N7G1mpyeP44jjwOfnDgG10rvE1/E4M1vbzHYALuR/e5GzgHnAzmbW18x6FXqgZGjjSuByMzvQzAYm75dDzew0+PowvauBN5M2jCX28q/Ke7j5wFVmtrGZbUIs9C8Tx6Nb4t/AMDPbLCm6E8jpgCTeATZO2l9rZoW+Sf4deBq43sy2MrP1knYsA1zawkzSIO1B9nK+UOCokgLLDCaOq35FPOJgON/cObkrsUe0iPY5HPAE4A7iH/FM4JS8Zb6xPZo/HPALGjkcMO+xt6aJQ+iSZaYSC1XDUTJzk+vL5SyzHfEog4XEIrJ3/vOWs+y/gH83sq2Dk7YuTp63s0gOB0zuH0Urd04C+ySv6cIkw3bk7fgFDiIWuCU0fThgF+BU4rDEYmLRfwT4SXL/aSRHlOSsMyh57n6a+3oABxB3/i4iFs61G2tvofYnt68B3Ed8384EzgQmAlNzlvlW8v6YR8sOB3yEwocD9svL0OxO+2q9NBwiJBXCzALxeNtr087SGDObSiz4h7fDY9UQe79/DCGc39bHK2fJT/MPCCEMTDuLlJZ2TkpZMrMMccz7SOLOuivSTSTScVS4pVz1Jw5BzAQOCSF8mXIekQ6joRIRkTKjo0pERMpMqYZK1I0XEWm5on4sV7Ix7hkzWjetc21tLbNmzWrnNJ2b2lwd1ObK15b21tXlHybfOA2ViIiUGRVuEZEyo8ItIlJmVLhFRMqMCreISJkp6qgS59x04mQ2S4F67/0mpQwlIiKNa8nhgNt776vnuB4RkU6q0wyVhGyW7N2eJW+9kXYUEZFOrai5Spxz7xAntw/AZd77/Anrcc6NJp7wFe/9xosXt2ye+uy8OXx24sFYJsPK511FpkfPFq1fzmpqaqivr087RodSm6tDtbW5Le3t1q0bFPnLyWILd533foZzbhXgAeA4733+SVhzhdb8cjK88ybZc0+HdTckc+wvsUyn+UJQUtX26zJQm6tFtbW5HX45WVThLqoyeu9nJP9+AtxOPF9cu7MBg+lx6PHw8rOEe28pxSZERMpes4XbObe8c65Hw/+BnYFXShVo2V2HY5ttR7jzesLrL5ZqMyIiZauYHveqwOPOuReJJ/2823v/t1IFMjPsoGOg7+pkLz+P8PlnpdqUiEhZavZwQO/928AGHZDla9Z9GTJHjSV71klkLzuHzMlnYTWFTiAtIlJ9Ou3eP1utH3bwcfDWG4RbJqUdR0Sk0+i0hRsgs+nW2A57ER76K9lnHk87johIp9CpCzeA7TMK1l6HMPnPhJkfpB1HRCR1nb9w13QlM/pU6NaN7KVnExYuSDuSiEiqOn3hBrCVa8kccTJ89CHhmr+gM9OLSDUri8INYOtugA0bQXj6EcLUe9OOIyKSmrIp3AC22z6w/iaEm64gvPNm2nFERFJRXoU7kyFz2Imw0spkx48jzJ2TdiQRkQ5XVoUbwJbvQeao02HOF2Qnnk/ILk07kohIhyq7wg1gaw7E9h8Nr/6LMMWnHUdEpEOVZeEGsG12wbbYnjDlRsIrz6UdR0Skw5Rv4TbDRh4Ndf3JXvFHwmefph1JRKRDlG3hBrDu3ckcNRayS8ledg5hyZK0I4mIlFxZF24AW7WOzKifwjtvEm6emHYcEZGSK/vCDWDf2xLb+ceEh+8h+9QjaccRESmpiijcADb8IBg0hHD1xYQP30s7johIyVRO4e7SJU5GtcyyZMefTVg4P+1IIiIlUTGFG8BWWpnM6FPg45mEyRdrMioRqUgVVbgB7NvrY8MPJDz7OOGhv6YdR0Sk3VVc4QawXYbDht8n3HIVYdrraccREWlXlVm4zcgccjys3IfsZecS5nyRdiQRkXZTkYUbwJZbgcyY0+GruWSv0GRUIlI5KrZwA1j/b2Ejx8DrLxLuvCHtOCIi7aKiCzdAZqsdsa13ItzjCS89k3YcEZE2q/jCDcQpYNcYQHbiBYRPP0o7johIm1RH4e6WTEZFIDv+HMKSxWlHEhFptaoo3ADWpy+ZQ0+E994i3Hh52nFERFqtago3gG2wGbbb3oRH7yP7z4fSjiMi0ipVVbgBbNgB8O31CddeSvjgnbTjiIi0WPUV7i5dyIw+GZZbgeyl4wjzv0o7kohIi1Rd4Qawnr3IHHkqzPqY7OSLNBmViJSVogu3c66Lc+5fzrkppQzUUWzQEGzvUfD8E4QH7kg7johI0VrS4z4eqKgZm2ynYfC9LQm3Tia8+UracUREilJU4XbO9QP2AK4obZyOZWbxfJV9ViM74Q+ELz9PO5KISLNqilzuT8CpQI/GFnDOjQZGA3jvqa2tbV2gmppWr9taS8aOY/aph9Plqj/R68wLsS7FPi3tI402p01trg7V1uaOaq81t2POObcnsLv3/mjn3FDgZO/9ns08bpgxY0arAtXW1jJr1qxWrdsW2ScfJky8ANtlOJl9RnXottNqc5rU5upQbW1uS3vr6uoArJhlixkq2Qr4oXNuOnAj8APn3LWtStaJZTbfHttuV8J9txFeeDLtOCIijWp2TMB7PxYYC5DT4z6gxLlSYfseQZg+jeyVF5L55ZrYKqulHUlE5Buq8jjuxljXrmTGnAZm8cc5ixelHUlE5BtaVLi991OLGN8ua1a7KpnDT4IPpxOuH592HBGRb1CPuwBbf2NsD0f4x0NkH7s/7TgiIv9DhbsRttd+MGRDwvWXEd57K+04IiJfU+FuhGW6xCGTHivG8e6v5qUdSUQEUOFukvVYMU5G9flnZK/6EyGbTTuSiIgKd3Ns7XWwnxwKLz5NuO+2tOOIiKhwF8N+sAe26TaE268lvPFS2nFEpMqpcBfBzLCDjoW+q8fJqD7/LO1IIlLFVLiLZMssS+ao02HxIrITziXU16cdSUSqlAp3C9hqa8Se97TXCbdNTjuOiFQpFe4Wymy2Lbb9HoQH7iQ894+044hIFVLhbgVzh8KAwWQnXUT46MO044hIlVHhbgWrSSajqulKdvw4wqKFaUcSkSqiwt1KtnIfMkecBDPeI1z7F50pXkQ6jAp3G9iQjbAf7k94cirhkb+lHUdEqoQKdxvZ7g7W25hw0+WE6f9JO46IVAEV7jayTIbMYSdCz15kx59DmDcn7UgiUuFUuNuBrdCTzJjT4cvZZCdeoMmoRKSkVLjbiQ0YhO17BLzyHOEen3YcEalgKtztyLbbFdt8KOGuGwiv/SvtOCJSoVS425GZYQccDautQfby8wmzP007kohUIBXudmbdl4mTUdUviTsr65ekHUlEKowKdwlY335kRv0U3nmTcPNVaccRkQqjwl0itvFW2I7DCH+fQvbpR9OOIyIVRIW7hGzvg2HguoSrLybMfD/tOCJSIVS4S8hqasiMPhW6dY9nil+4IO1IIlIBVLhLzHr1JjP6FPjow9jz1mRUItJGKtwdwNb5LvajkYRnHiM8fHfacUSkzKlwdxDbdW/YYDOCv5Lw1htpxxGRMqbC3UEskyFzyAnQqzfZy84lzP0y7UgiUqZUuDuQLb9C/HHO3C/JXnE+Ibs07UgiUoZqmlvAObcM8CjQPVn+Fu/9r0sdrFJZ/7WxEUfGHZV/vQkbNiLtSCJSZorpcS8CfuC93wDYENjVObd5aWNVNtt6J2yrHQhTbiS8/FzacUSkzDTb4/beB2BecrVrctExbW1gZjBiDOHdt8lO/CNLv/NdyHRNO5aIlAkr5rhi51wX4DlgIHCJ9/60AsuMBkYDeO83Xrx4casC1dTUUF9f36p1y039zA+YffIh1NT1p9fvL8W6dks7Uoeppte5gdpc+drS3m7dugFYMcsWVbgbOOdWAm4HjvPev9LEomHGjBlFP26u2tpaZs2a1ap1y1F44Umyl/weG7obmZFHpR2nw1Tb6wxqczVoS3vr6uqgyMLdoqNKvPdfAFOBXVucSgqyDTdnuR+PJEy9l+yTD6cdR0TKQLOF2znXJ+lp45xbFtgR0C9I2tEKI4+EwesRrrmE8OG7accRkU6umB73asDDzrmXgGeAB7z3U0obq7pYl5o4n8myy8fJqBbMTzuSiHRixRxV8hKwUQdkqWq2Yi8yo08he/4vyU66iMyY0+LRJyIiefTLyU7EBq+HDT8Ynv8n4cG70o4jIp2UCncnYzv/CDbanHDrJMK019KOIyKdkAp3J2NmZEYdD71XiZNRzfk87Ugi0smocHdCttzycTKq+fPITjiPsFSTUYnIf6lwd1LWbwA28ij498uEO69LO46IdCIq3J1YZssdsG12Jtx7C+HFp9OOIyKdhAp3J2f7j4b+a5OdeAHh04/SjiMinYAKdydnXbuRGXMamJEdP46wpHWTd4lI5VDhLgPWpy+Zw06E994m3DAh7TgikjIV7jJh390U2/0nhMfuJ/uPB9OOIyIpUuEuIzZsBKy7AeG68YT33k47joikRIW7jFimC5nDT4Lle8Tx7vnzml9JRCqOCneZsZ4rkTnyVJj9KdmrLqQlJ8IQkcqgwl2GbOC62D6HwAtPEe67Le04ItLBVLjLlO2wF7bJ1oTbriH8u6mzyIlIpVHhLlNmhh18LKy6GtkJ5xK+mJ12JBHpICrcZcyWWY7MmLGwcEEs3lV0Nm2RaqbCXeZs9f7YgcfAf14j3H5N2nFEpAOocFeAzOZDsaG7E+6/nfD8E2nHEZESU+GuEOYOgwGDyU66kPDxjLTjiEgJqXBXCOvalcyRp0GXLmQvPZuwaFHakUSkRFS4K4j17kPmsJNgxnuE6y7Vj3NEKpQKd4Wx9b6H7bkv4Ym/Ex67P+04IlICKtwVyPbcF4ZsRLhhAuHdaWnHEZF2psJdgb6ejKrnimQvHUf4am7akUSkHalwVyjr0ZPMmNPhi9nxtGfZbNqRRKSdqHBXMBswGNv3MHj5WcK9t6QdR0TaiQp3hbOhu2ObbUe483rC6y+mHUdE2oEKd4UzM+ygY6Dv6mQvP4/w+WdpRxKRNlLhrgLWfRkyR42FxYvJXnaOJqMSKXMq3FXCVuuHHXwcvPUG4dZJaccRkTaoaW4B59wawNVAXyALTPDeX1jqYNL+MptuTfat1wkP3kX2W+uQ2XTrtCOJSCsU0+OuB07y3q8LbA4c45wbUtpYUiq2zyhYex3C5D8TZn6QdhwRaYVmC7f3fqb3/vnk/3OB14HVSx1MSsNqupIZfSp07RrPFL9oYdqRRKSFrCUTETnn1gIeBdbz3s/Ju280MBrAe7/x4sWLWxWopqaG+irbeZZGmxe99CxfnHkCy2y9Iz1P+DVm1qHb1+tcHaqtzW1pb7du3QCK+kMsunA751YAHgHO8t43d2rxMGNG6+aErq2tZdasWa1at1yl1ebs3Z5wx7XYiDFktt+9Q7et17k6VFub29Leuro6KLJwF3VUiXOuK3ArcF0RRVvKhO22D6y/CeGmKwjvvJl2HBEpUrOF2zlnwETgde/9H0sfSTqKZTJkDjsRVlqZ7PhzCHPnNL+SiKSumB73VsCBwA+ccy8kl479Xi0lY8v3IHPU6TDnc7ITzydkl6YdSUSa0exx3N77xyly3EXKk605ENt/NOGavxCmeOyH+6cdSUSaoF9OCgC2zS7YFtsTptxIeOX5tOOISBNUuAVIJqMaeTTU9Y9DJp99mnYkEWmECrd8zbp3j5NRLV0aJ6NasiTtSCJSgAq3/A9btY7MqJ/CO28Sbp6YdhwRKUCFW77BvrcltvOPCA/fQ/apR9KOIyJ5VLilIPvxQTBoCOHqiwkfvpd2HBHJocItBVlNTZyMapllyY4/m7BwftqRRCShwi2NspVWJjP6FPh4JmHyxbRkQjIRKR0VbmmSfXt97McHEp59nPD3KWnHERFUuKUItutw2PD7hJuvJEx7Pe04IlVPhVuaZWZkDjkeVu5D9rJzCXO+SDuSSFVT4Zai2HIrkBlzOnw1l+wVmoxKJE0q3FI06/8tbOQYeP1Fwp03pB1HpGqpcEuLZLbaEdt6J8I9nvDSM2nHEalKKtzSYrb/aFhjANmJFxBmfZx2HJGqo8ItLWbdksmoQohnzlnSuhNDi0jrqHBLq1ifvmQOPQHenUa48fK044hUFRVuaTXb8PvYbnsTHr2P7D//nnYckaqhwi1tYsMOgG+vT7juL4QPpqcdR6QqqHBLm1iXLmRGnwzLrkD20nGE+V+lHUmk4qlwS5tZz15kjjwVZn1EdvJFmoxKpMRUuKVd2KAh2N6j4PknCA/ckXYckYqmwi3txnYaBt/bknDrZMKbr6YdR6RiqXBLuzGzeL7K2r5kJ5xL+PLztCOJVCQVbmlXtuxyZI46HRZ8RXbCHwhLNRmVSHtT4ZZ2Z/3Wwg44Bt58hXDHtWnHEak4KtxSEpkttse225Xwt1sJLzyZdhyRiqLCLSVj+x4Baw4ke+WFhE9mph1HpGKocEvJWNeuZMacBmbxxzmLF6UdSaQiqHBLSVntqmQO/xl88A7h+svSjiNSEWqaW8A5dyWwJ/CJ93690keSSmPrb4LtuS9hyk1k116HzDY7px1JpKwV0+OeBOxa4hxS4Wyv/WDIhoTrLyO891bacUTKWrOF23v/KDC7A7JIBbNMFzKHnwQ9ViQ7/hyyX81NO5JI2Wp2qKRYzrnRwGgA7z21tbWtC1RT0+p1y1XVtLm2lsWn/Z7Pf3k0cy/6Hb1POxvLVM9ulqp5nXNUW5s7qr3tVri99xOACcnVMGvWrFY9Tm1tLa1dt1xVVZt798X2OZSFN05g0XUTyOy2T9qJOkxVvc6JamtzW9pbV1dX9LLV092RTsN+sAfdt96RcPu1hDdeSjuOSNlR4ZYOZ2b0PPp0WLUuzmfyxWdpRxIpK80WbufcDcATwLedcx845w4rfSypdJmGyagWLyJ72bmE+vq0I4mUjWbHuL33+3dEEKk+VtcfO+hYwuXnEW6bjKlPIFIUDZVIqjKbbYttvwfhgTsJz/0z7TgiZUGFW1Jn7lAYMJjspAsJH32YdhyRTk+FW1JnNclkVDU1ZMePIyxamHYkkU5NhVs6BVu5D5kjToYZ7xGuvVRnihdpggq3dBo2ZCNsr/0JTz5MePS+tOOIdFoq3NKp2B4O1tuYcOMEwvT/pB1HpFNS4ZZOxTIZMoedCD17kR1/DmHenLQjiXQ6KtzS6dgKPcmMOR2+nE124gWEbDbtSCKdigq3dEo2YBC27+HwynOEe25OO45Ip6LCLZ2WbbcbtvlQwl3XE177V9pxRDoNFW7ptMwMO+BoWG0NspefT5j9adqRRDoFFW7p1Kz7MnEyqvolyWRUS9KOJJI6FW7p9KxvPzKjfgpv/5twy6S044ikToVbyoJtvBW24zDCQ38l+/SjaccRSZUKt5QN2/tgGLgu4eqLCTPfTzuOSGpUuKVsWE0NmdGnQrfuZC8dR1i4IO1IIqlQ4ZayYr16x8moPvow9rw1GZVUIRVuKTu27gbYj0YSnnmM8PDdaccR6XAq3FKWbNe94bubEvyVhLfeSDuOSIdS4ZayZJkMmUNPhF694/Hdc79MO5JIh1HhlrJly68Qf5wz90uyV5xPyC5NO5JIh1DhlrJm/dfGRhwJr71A+OtNaccR6RAq3FL2bOudsC13INx9E+Hl59KOI1JyKtxS9swMGzEGVl+L7MQ/Ej77JO1IIiWlwi0Vwbp3J3PUaZBdGn+cs0STUUnlUuGWimGr1JE55AR4dxrBX5F2HJGSUeGWimIbbY7tMpww9V6yTz6cdhyRklDhlopjPz4QBq9HuOYSwofvph1HpN2pcEvFsS5dyIw+BZZdPo53L5ifdiSRdqXCLRXJVuwVi/enM8lOvkiTUUlFUeGWimWD18OGHwzP/ZPw4F1pxxFpNzXFLOSc2xW4EOgCXOG9H1fSVCLtxHb+EeGt1wm3TiIMGIQNHJJ2JJE2a7bH7ZzrAlwC7AYMAfZ3zundL2XBzMiMOh56rxIno5rzRdqRRNqsmB73ZsA07/3bAM65G4FhwGulDCbSXmy55ckcdTrZs08he+ZPYfkeqeSYVVPD0vr6VLadlmpr8+xeK8OJvy35doop3KsDuSf4+wD4fv5CzrnRwGgA7z21tbWtC1RT0+p1y5Xa3AFqa1l02jgWPDQFUtpRaRmjJltdO0mrrc1devSkRwe8r4sp3Fbgtm+8Et77CcCEhvtnzZrVqkC1tbW0dt1ypTZ3kP4D4ZATOnabOfQ6V76V2tDeurq6opct5qiSD4A1cq73A2a0MJOIiLSTYnrczwCDnHMDgA+B/YARJU0lIiKNarbH7b2vB44F7gNejzf5V0sdTERECivqOG7v/T3APSXOIiIiRdAvJ0VEyowKt4hImVHhFhEpMyrcIiJlxko03WX1/FRKRKT9FPrB4zeUqsdtrb04556wRktoAAAKpElEQVRry/rleFGbq+OiNlf+pR3aWxQNlYiIlBkVbhGRMtMZC/eE5hepOGpzdVCbK1+HtLdUOydFRKREOmOPW0REmqDCLSJSZoqaZKqUknNaPgt86L3fM5k+9kZgZeB54EDv/eI0M7a3Am2+DtgEWAI8DRzpvV+SZsb2lt/mnNv/DBzivV8htXAlUuB1NuB3wE+ApcCl3vuL0szY3gq0eQfgD8RO4jxglPd+WpoZ25Nzbjowl/h61nvvN3HOrQzcBKwFTAec9/7z9txuZ+hxH0+cLrbBOcAF3vtBwOfAYamkKq38Nl8HrAOsDywLHJ5GqBLLbzPOuU2AldKJ0yHy2zyKeFKSdbz36xI7KJUmv82XAiO99xsC1wO/TCVVaW3vvd/Qe79Jcv104KGkhj2UXG9XqRZu51w/YA/giuS6AT8AbkkWmQz8KJ10pZHfZojT5nrvg/c+EHvc/dLKVwqF2pz0zP4AnJpWrlIq1GbgKOA33vssgPf+kzSylUojbQ5Az+T/K1IdZ88aRqxdUKIalnaP+0/EP9xscr038EVy8gaIp01bPY1gJZTf5q8557oCBwJ/6+hQJVaozccCd3nvZ6YTqeQKtXltYF/n3LPOuXudc4PSiVYyhdp8OHCPc+4D4nt7XBrBSigA9zvnnktOmA6wasP7Ovl3lfbeaGqF2zm3J/CJ9/65nJsL/eSzYo5XbKTNuf4CPOq9f6wDY5VUoTY75+qI47x/Ti1YCTXxOncHFiZfqS8HruzwcCXSRJtPBHb33vcDrgL+2OHhSmsr7/33gN2AY5xz23bERtPscW8F/DAZ3L+ROETyJ2Al51zDTtNKOzHxN9rsnLsWwDn3a6AP8LP04pVEodf5VWAgMC25fTnnXMXssKLx1/kD4NZkmduB76YTryQKtfluYAPv/VPJMjcBW6aUryS89zOSfz8hvqabAR8751YDSP5t9yGx1I4q8d6PBcYCOOeGAid770c6524G9iG++AcDd6aVsb010uYDnHOHA7sAOzSMf1aKRtq8Z+4yzrl53vuBKcQriSZe53HED64rge2AN1ML2c4KtZk4tvuRc26w9/5NYCfydlCXM+fc8kDGez83+f/OwG+Au4i1axwlqmFpj3EXchrws6QH1huYmHKejjAeWBV4wjn3gnPuV2kHkpIYB+ztnHsZOJvKPHroa8m+qiOAW51zLxLHuE9JN1W7WhV4PGnb08Dd3vu/EV/nnZxz/yF+WLX7uL5+8i4iUmY6Y49bRESaoMItIlJmVLhFRMqMCreISJlR4RYRKTMq3FKQc26Uc+7xNqw/yTn3uyKW6++cm5fMXVJ2nHNTk+PwS/HY2zjn/l2Kx5bylvq0rtI+nHMBGJQ7ZaZz7gxgoPf+gNSCNcN7/x7QoVO6Jj8QuTb5GXanlUx98O1SbsM51414TPm+xJkaPwXu8N6fWMrtStuoxy3trlx7z51JzrQPpTaWOBf8ZkAPYHvgXx20bWkl9birREMvE7iA+OvUpcDPvfdXJff3Jk4CNBR4A7gvb/11iJNCbUzslf2f994n900CFgBrEn/KPSxv3VeAsd77vybXuwIzgR2BL4B3gK7e+3rn3FTgMeJPw78LPAGM8N7PStY9CPgtsZf+J+J87Yd77x8s0ObdgfOIc2DPSdp+KXAv0N05Ny9ZdDDwGXEueJfc5oHTvPeLkscaBpwJfCtp/zHJr+Ryt7da8rxd7b0/r0Ce6cBlxF8QrgbcARzlvV+Y8/r8mTgx0wPOuYnkfDNwzq0BXAhsQ+x03eC9Pza571DirxL7En/FN9p7/25+hgI2BW5vmHODOPH/9CLWkxSpx11d+hLnRF6dWPAucc71Su67BFhILCiHJhfg6zkZHiBOhL8KsD/wF+fcd3IeewRwFrHXlj82fjWQO1yzOzDTe/9CIzlHAIck2+pGnPcC59wQ4gyKI5OcDW1pzETi2YR6AOsBf/fef0WcyW2G936F5DID+AWwObAhsAGxB/rLZLubJW04hTicsC15xc05txbwCHBxoaKdYyRxXpq1iR8YuScW6Es889OawOjclZJvMVOAd4lnVlmd5EQMzrkfAT8HhhMnKnsMuKGJDLmeJE4xcbRzbv1kTnzp5NTjri5LiBP51xPnSJ4HfNs59wywN7B+Uthecc5NJhYogD2B6Q29c+B559ytxMnAXk1uu9N7/4/k/wuda+i4ArEn+X/OuZ7e+znEHuc1TeS8KpmUCOecB36Y3L4P8Ffv/ePJfb8CftpMe4c4515MTh3V1OmjRgLHNZzcwDl3JrF3/H/ED7krvfcPJMt+mLfuEGIBHuu9b65gXuy9fz/ZxlnEHnZD8c4Cv87p5eeutxlQB5ySM199wwfkkcDZ3vvXk/V+D/zcObdmEb3us4nPy0jiN5LPnHNjvfeTm15N0qQed+VYCnTNu60rsXg1+Cznjx5gPnHIoQ/xQ/z9nPty/+DXBL7vnPui4UL8Q++bs0zuuv8j6dH+gzjB0krEHu91TbTlowIZIRaur7fjvZ9PHOJozN7E3v27zrlHnHNbNLFsHf/b5neT2yAOtbzVxLojicX8liaWaZD/HNflXP/Ue7+wkfXWAN7Ne/0arAlcmPPazCbObd/sSUi890u995d477cifps4C7jSObduEW2RlKjHXTneI36Fzp02cwDFTR36KVBPLA5vJLf1z7n/feAR7/1OTTxGc7OVTSbOhlcDPOG9z++1FmMmOUdZOOeWJc4gWZD3/hlgWDKmfixx3HqNRrLOIBbAhm8Q/fnvXPDvE4c2GnMGsCtwvXNuP+/90iaWXSPn/7nboJFcDd4H+jvnagoU7/eBs7z3TX0YNst7v4A4fHYm8VtExUzBWmlUuCvHTcAvkylDZxB37u0FNNXLBGKvyzl3G3BGspNrLeI8wtOTRaYA45xzB/LfE9xuCMxr+HpehDuI49OrAucWuU6+W4AnnXNbEs8kfiaFz5rUcJjbT4Ap3vsvnXNziN9KAD4GejvnVvTef5ncdgPx+XuGWEB/RRzigThWfr9zbgrwMHF8vYf3vuFDbkmyrTuAa5xzBzQxr/oxyePMJ45L31Rk258mfnCNS066sRTYOBmeGg/81jn3gvf+VefcisDO3vubk+diKjDVe39GgefpBOAF4KmkHSOJ+yl0ZEknpqGSyvEb4J/Ecc/PicVxpPf+lSLXP5Y4JPERMIl4hAkA3vu5xEni9yN+KHxEPAKje7Hhkt7crcRvAbcVu17eY7wKHEf88JgJzCWeXWRRI6scCExPivYYkh2kScG9AXg7GV6oA35H/DB4CXgZeD65De/908SdpRcAXxJ3Qq6Zl20xcefgKsShhsb+tq4H7gfeTi7N/kgpefylxA/igcRvVx8Qj73Ge3878fW4MWnrK8ThqAZrEIeqClkAnE98TWcBxwB7e+/fLiaXpEPzcUuHSXYmDm6vHwQ551YgHk44yHv/Tns8ZiklhwMWPHSxhNvsB9zsvW/2m5eUDw2VSIdwzq1MPDrjwDY+zl7AQ8QhkvOIvePpbc1Xqbz3H1DEcJmUFw2VSMk5544g7kC713v/aBsfbhhxuGYGMAjYz3uvr41SVTRUIiJSZtTjFhEpMyrcIiJlRoVbRKTMqHCLiJQZFW4RkTLz/73cRjJZhDLoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7a6fe9887550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use a lambda for a call payoff function\n",
    "# equivelant to:\n",
    "#\n",
    "# def call_payoff(S, K):\n",
    "#     return np.maximum(S - K, 0.0)\n",
    "call_payoff = lambda S, K: np.maximum(S_ - K, 0.0)\n",
    "\n",
    "# and put payoff function\n",
    "put_payoff = lambda S, K: np.maximum(K - S_, 0.0)\n",
    "\n",
    "# plot the call payoff \n",
    "plt.figure(1)\n",
    "plt.title('Call option payoff at expiration')\n",
    "plt.xlabel('Underlying stock price, S')\n",
    "plt.plot(S_, call_payoff(S_, K))\n",
    "\n",
    "# plot the put payoff\n",
    "plt.figure(2)\n",
    "plt.title('Put option payoff at expiration')\n",
    "plt.xlabel('Underlying stock price, S')\n",
    "plt.plot(S_, put_payoff(S_, K))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# (Some of) the math\n",
    "\n",
    "The famous Black-Scholes model:\n",
    "\n",
    "> $$\\frac{\\partial V}{\\partial t} + \\frac{1}{2}\\sigma^{2}S^{2}\\frac{\\partial^{2} V}{\\partial S^{2}}+rS\\frac{\\partial V}{\\partial S}-rV=0$$\n",
    "\n",
    "The key financial insight behind the equation is that one can perfectly hedge the option by buying and selling the underlying asset in just the right way and consequently create an arbitrage free environment. This hedge, in turn, implies that there is only one right price for the option, as returned by the Black–Scholes formula. For this dynamic hedging to work, one must continuously hedge the movement in the underlying stock with the option which in reality, is prohibitive due to transaction costs.\n",
    "\n",
    "I'll skip the derivation of the option pricing formulas from the partial differential equation but the result solves for a call option\n",
    "\n",
    "> $$C(S, t) = N(d_{1})S - N(d_{2})Ke^{-r(T-t)}$$\n",
    "\n",
    "and put option\n",
    "\n",
    "> $$P(S, t) = N(-d_{2})Ke^{-r(T-t)} - N(-d_{1})S$$\n",
    "\n",
    "where\n",
    "\n",
    "> $$d_{1} = \\frac{ln(\\frac{S}{K})+(r+\\frac{\\sigma^{2}}{2})(T-t)}{\\sigma\\sqrt{T-t}}$$\n",
    "\n",
    "and\n",
    "\n",
    "> $$d_{2} = d_{1} - \\sigma\\sqrt{T-t}$$\n",
    "\n",
    "and finally (as we discussed above)\n",
    "\n",
    "* $N(\\cdot)$ is the cumulative distribution function of the standard normal distribution\n",
    "* $T - t$ is the time to maturity (I use $t$ in the code)\n",
    "* $S$ is the spot price of the underlying asset\n",
    "* $K$ is the strike price\n",
    "* $r$ is the risk free rate (annual rate, expressed in terms of continuous compounding)\n",
    "* $\\sigma$ is the volatility of returns of the underlying asset (I use $vol$ in the code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Black-Scholes formula in Python\n",
    "\n",
    "Finally, enough math and some code. We'll be using Numpy so we vectorize the function which will help us later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N(z):\n",
    "    \"\"\" Normal cumulative density function\n",
    "\n",
    "    :param z: point at which cumulative density is calculated \n",
    "    :return: cumulative density under normal curve\n",
    "    \"\"\"\n",
    "    return norm.cdf(z)\n",
    "\n",
    "\n",
    "def black_scholes_call_value(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes call option\n",
    "\n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: BS call option value\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    d2 = d1 - (vol * np.sqrt(t))\n",
    "    \n",
    "    return N(d1) * S - N(d2) * K * np.exp(-r * t)\n",
    "\n",
    "def black_scholes_put_value(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes put option\n",
    "\n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: BS put option value\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    d2 = d1 - (vol * np.sqrt(t))\n",
    "    \n",
    "    return N(-d2) * K * np.exp(-r * t) - N(-d1) * S\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Black-Scholes call value 3.20\n",
      "Black-Scholes put value 2.79\n"
     ]
    }
   ],
   "source": [
    "print ('Black-Scholes call value %0.2f' % black_scholes_call_value(S, K, r, t, vol))\n",
    "print ('Black-Scholes put value %0.2f' % black_scholes_put_value(S, K, r, t, vol))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we vectorized the function, we can plot the Black-Scholes call value along a range of underlying stock prices. Going back to our example of the option value at expiration, let's see what it looks like with six, three and one month left to expiration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the value of the option with six months to expiration\n",
    "black_scholes_call_value_six_months = black_scholes_call_value(S_, K, r, 0.5, vol)\n",
    "\n",
    "# get the value of the option with three months to expiration\n",
    "black_scholes_call_value_three_months = black_scholes_call_value(S_, K, r, 0.25, vol)\n",
    "\n",
    "# get the value of the option with one month to expiration\n",
    "black_scholes_call_value_one_month = black_scholes_call_value(S_, K, r, 1.0/12.0, vol)\n",
    "\n",
    "# get payoff value at expiration\n",
    "call_payoff_at_expiration = call_payoff(S_, K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7a6fe1f5cf28>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEaCAYAAAAxEb1YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXd8HMX5/99z6tW2ii25yRUX2bj33isuMlxCCQESjL8QvrQEEpJ8Q8IP0hMIkEJCrzmK3DuuyF22ZWRbLtiWi3qvd6e7nd8fezJCVjnJOt2dNO/XSy/79mZ3PrO79+zsM888I6SUKBQKhcJ7MLhbgEKhUCiahjLcCoVC4WUow61QKBRehjLcCoVC4WUow61QKBRehjLcCoVC4WV4heEWQuwSQvynBY/XSwghhRCTW+qYTaj7OSHE+Vaq6z4hhK016moOQojpjuvQ3QO09BBCfCGEKBdCuCVGVghxSQjxC3fU7ai/RX9nrqIldXrSPdgU3Gq4hRBvO05a9V+xEGK/EGKhO3U1hBDifiFEshCiRAhRKoQ4LYT4t7t1eSn7gFggw91CgGeBzsBwdE0uQwjxHyHErjq+GgP81ZV1O+q/x10PJ3cihLAJIe6rtdmT7kGn8YQe9170ExcLjAeOAquFEH3dqqoOHBf9n8BbwGhgFPBTwMeNsrwSIYS/lNIqpcySUmru1gP0Bw5JKc9JKbPcIUBKmSulLHdH3c1FCOEnhBDu1tFcPOwedB4ppdv+gLeB7bW2hQESWF5j2y7gPzU+z3FsKwCKgd3A2FrHCQVeAq4AFuAS8Kzju16OOibXKP8MUAhMb0DvauBTJ9o1CtgMlABlwCFgnOO754DzwFIgDSgHdgJ9a+x/H2ADJqE/yCqAw8CoWvWMB/YAlQ7tHwKdax+nDm1bHbpygc+BuBrfdwc+A/Icx70A/KSBtk53nMvbHO00AyeBOXWUWQR86Sjzoxrbu9co2xf4xHFtK4ATwGJn9dejMQz4l6O8GTgCzK3xvaz193YDx1oIJDvuqRzg70BI7XsaeBK45mjDZ0BUjetfu777HN9dAn7RBN29HPsbgXWOui4A33Piet3QXhy/M+CXQJbjGrxdT/sedejV0H9rfsDvHG22AqeAu2rVLYF7am3bXvN8A5GO618OZAPPA+9Qw044o7OOdl+q3e5a56N7rc8Lgf3ov4FkIN7x96XjPB8CBjflt9WSf57Q476OEMIfeBD9R3G0gaKhwGvohmsicA7YLISIdBxHAOuBJeg32CDgXvSTWbtOgxDiFeAxYJqUclcD9WYCo4UQtzTQhnh0Y1oIzARGoL/+1jzXscD/AHc79HcE3qx1KAPwW4eukY7jmYQQvo56YtBvkqvAWHTDOQTdSNSnbTD6Q24/+hvDTMAObBNCBDqK/R3oAMxGP28/cNTRGH8BfuNo7wFgrRCiW60yfwb+4Dju6jr0xaC/unZCv3ZD0X+cWhP018WbwDzgHoe+JGC9EGKg4/tYxzE/dPz/sboOIoS4FViLfn2HA98HFqO/hdVkLLoBmI9uAG7lm+v7J0c9+/nmTfO/zdRdze+A9xz1mIC3hBD96znmPvSHZnW7a7f3diDCof8uYBnwdB3tm+n4bhj6Q+VF9N/u4+j34fvA+0KIWfXoqI+3HMdc7Kiju6Oe2jijsyZj0O+Vx/mm3Q3xAvBzdGNsBT4C/gH8qsa2t6oL38S92Txc8TRw9g/9KWlDf0KVof9AywBjrXK7qNHjruM4BnTDdrfj8yz0p+boesr3cnw/C/3pnoYTT0YgBv2JK9Gf4P8FVgLBNcq8B6QAhnqO8ZyjzdE1tn3X0fZAx+f7HHWMrFFmvGPbAMfn59ENqn+NMsMcZabWOI6txvdvAx/X0hOA3oNY5vicAjzXhGs43VHnD2ps8wXSgf9Xq8z36tm3e402ZVFPz8kZ/XXs089Rx8Ja248Cbzp7j9W4todqbVvquHZxNTSWAR1qlJnr0NDf8fk/wK46jn8JR4/bGd18cx8/WevclwEPNdCOe3D0OOv4nZ2ote2fwP5a16AICK2xLRi9s/VwrX0TgR01PjfY40Z3V0lgVo3v/dDfmmv3uBvUWU+7bTjebhq4B6s/L6tR5g7HthU1ti13bAtt7r15M3+e0OM+iN57GY7es/wN8I4QYl59Owghegsh3hNCnBdClKC7JDoAcY4io4BCKeWRRup+C71XN0lKmV7j+FOEEGU1/p4FkLovbDIwGL03XI7egzwphOhco+4vZMM+swwpZc3e/zVAoA+OVSPRjWjNMgBdHP/GAweklNbrO0iZgu46iq+n3jHA8pptA/KBQPQfDejupWeFEAeFEL8XQkxtoB012V9Dhw3Hq2StMocaOcYoYJ+s38/rjP7aVGvYU2v7Huo/T/VR/TZVk93o165mW09JKYtrfE5y/DuoCXU1Rffx6v84zn0239wnTeV4rc/X6jjWaSllWY3P/QD/OrTurkNrQ1S3+UD1BillFbqLqDk6b4aav73qMY8TdWyr/s02595sNr4tfcBmUCmlrBked9zxevVzYEs9+6xH98E+gv40tqL3hP1rlJFO1L0BeAD9lfaDGtuPoD9IqimouZOU8jRwGviXEOJ54Cy66+PXTtZtrfW5unzNB6kmpbQ3Uqa+eurbbkDvNf6uju/yAaSUbwkhNqOfkxnAJiFEopTynnqOWR91DVg5M/DW0LlrVH8TEI3UVR9NPectTV2667qfmtspc+ZY9V3H2rpqa5XceF/4OXGcumjJNtdFVa1j17fNUOPflro3G8UTDHdd2NBfv27A4ccejP4KucWxrTvf7q0mAxFCiNGN9Lo/QO8VvCOE8JVSvgMgpaxEH0B0hkvor0PV9ScDs4UQhkZ63TfLSeD+6ugMACHEMPQ3j5P17HME3Q/6tXS8y9WFlDIT/W3kLSHERuAjIcTDUsqSBvSMRx+QwuGHH4Pu52wKycCDQoiQenrdTumvRfW5mApsrLF9CnCsifpOAtNqbZuG/iM+VWPbICFEeI3zNdHx72nHv1Yaj0RqSd21qb5ffGp1DprLeXRXyTS+fe9NrfU5B+ha/UEIEYD+W77o2FR9DicAXzjK+KK/iZ1tAZ3OnPfm0px7s9l4gqvEXwgR4/jrK4R4GH1AJrGe8oXog4wPCiFuEUJMQB84qKxRZgd6mOF/hRBLHa6VSUKIH9Y+mJTyY/TBjX8JIR5sSKgQ4h9CiF85XClxQohR6CPe4Xwz2PYH9FejD4QQox1tusOhsyV51VHv20KIIUKfTPQe8KWUcm89+7yI/rr+vhBirOO8zBBCvCyE6ONo46tCiIUO3fFAAvpbTWkjen7q2G8Q+iBOF8e/TeHv6PfkGsf16i2EWCyEWOCs/tpIKb9GH8f4uxBinhBioBDiZfQBtD82Ud8fgZFCiL84jjMfeAX4QEp5uWa1wLuO6zIVfSB9g5TynOP7i8BAIUS8ECLKYcBcqbs21YZyiRAiWggRejMHk1JWAH8Dnnfc6/0d7sWl6Nesmu3AKiHEBCHEEHS/sH+N45xDj455TQgxzTHg9y/0+7wljOFFYIYQoqsQIqoFjleTJt+bN4MnGO4p6NEamcBX6O6Pn6L7kG/A0Yu9Az1s7AT6xX/JsX91GYkeerYRfdDiDHrvr86LJaX8DD2k6m+OB0d9bEN/+n+E3gPYiD46vVBKuc1xrK/QBzii0Xvzx4Efo48wtxhSymz0Qa/u6KGC64FUYEUD+5xG7/2ForuhTgH/BoLQB5xAf5V9yXGsPUAIsMCJXsSP0QcXj6OHMS6VUjoTjVJTXyYwGf0hsRG9t/aCQ5Oz+uvih47y76P7LiehhximNVHfCfRol2mO47yH7m5bVavoIXTX3TZHvSeB+2t8/wb6NduH3gm505W662jHYeBl9N9GNnon4Gb5Ofq1eAm9vfegD0R+UaPMj9Hvqy3AJvT763Ct49zvKLMJfRDyGvp5NLeAxqfQf78XqSPC7Ga4iXuzWYhW6NUr2jBCiOnoceg9mmqo2yJCiLfRIxRmu1tLW0AI4YMe9bVWSvmUu/V4Cp7q41YoFO0Qh2upM7ofPwx4Aj3s8W33qfI8lOFWKBSehA/wC/QQwyp0t8kMhwtS4UC5ShQKhcLLaLTHbTQaB/DtKbl9gP8zmUwvuUyVQqFQKOqlST1uo9Hogz7KO85kMqU3UFR14xUKhaLpOJVpsak+7lnA140YbQAyMpqX3jYqKoq8vLxm7eutqDa3D1Sb2z43096uXbs2XshBUw33d9FjmG/AaDSuRE+4hMlkIiqqefHtvr6+zd7XW1Ftbh+oNrd9Wqu9TrtKjEajP/oqEfEmkym7keJS9bidR7W5faDa3PZpgR63U66SpsycXAAcdcJoKxQKhcKFNMVVcif1uEmcQUqJ2WxG0zREAysdZWdnY7FYmluNV1KzzVJKDAYDgYGBDZ4nhULRfnHKcBuNxmD05cIeam5FZrMZPz8/fH0brtLX1xcfn/a1hGPtNttsNsxmM0FBQW5UpVAoPBWnDLfJZKpAXwuu2Wia1qjRVuj4+vq2u7cOhULhPK2WHVC99jcNdb4UCkV9eEJaV4VCofB65PnTlCd+0HjBFqDdGO7i4mLefvttp8tfvnyZxYsXM2nSJFatWoXVWnulJLhy5Qp9+/Zlzpw5zJkzh2eeeaYFFSsUCm9Aahrapk/R/vgzKrckIs2Vje90k7Qbw11SUsK7777rdPkXXniBBx98kKSkJDp06MBHH9UdUBMXF8e2bdvYtm0bv//971tKrkKh8AJkSSHay79Gfv4uYsQEIv78NiLQ9UEF7cZwv/jii6SnpzNnzhyef/75BstKKUlKSmLRokUA3HHHHWzZUt+6xQqFoj0iT6eg/eZxOHcScc/DiIeexhByU6vAOY1bwjy0j/+NvHKx7u+EoDmpZkWP3hi+W/+Skc8++yxnzpxh27ZtlJWVMWfOnDrLvfbaa0RFRdGhQ4frUTCxsbFkZWXVWf7y5cvMnTuXsLAwnn76acaNG9dk7QqFwnuQdjty3UfIjZ9Al24YHv81onuvVtXQLuPzQkND2bZtW73f5+fn37CtriiPzp07c+jQISIiIjhx4gQPPPAAO3fuJCwsrEX1KhQKz0AW5KL9589w7hRi4izEXQ8hAgJbXYdbDHdDPWNfX19sNptL6y8rK2P58uV1fvfaa6/Rv39/iouLsdls+Pr6kpmZSZcuXW4oGxAQQECAvkD3rbfeSq9evbhw4QLDhg1zqX6FQtH6yKP70d55Bew2xA+ewDB+htu0tJsed0hICGVlZUDjPW6AiRMnsmHDBpYuXconn3zC3LlzbyiTn59Px44d8fHxIT09nYsXL9KzZ0+X6FcoFO5BWi3IT95E7toEcf0wPPhjRBfnU7C6gnYzOBkREcGYMWOYOXNmo4OTAD//+c95/fXXmTRpEoWFhdx5550AbN26lT/+8Y8AHDhwgNmzZzN79mxWrlzJb3/7Wzp16uTSdigUitZDXktHe/HHyF2bEHOXY/jp791utMF1a07ekNa1oqKC4ODgRndsDVeJp1FXm509X95Ke0v3CarN3oSUErl7M9L0BgQGYXjgCcSQkY3u11ppXduNq0ShUCicQZaXor37KhzdD/EjMDzwOCLcs96kleFWKBQKB/JsKtp//gIlRYg77kfMXooweJ5HWRluhULR7pF2O3LDf5HrTRDdRfdl9+rvbln1ogy3QqFo18j8XLQ3HLHZE2bosdmBnj2+pAy3QqFot8jkfWjvvgJ2ze2x2U1BGW6FQtHukBYL0vQGcs9m6NVfj83uHOtuWU7jeV53F+GKtK579uxh/vz5zJo1i/nz5/Pll19e/+72229nypQp11O+emNIlELRFtFjs59C7tmMmJeA4ZnfeZXRhnZkuF2R1jUiIoK3336bL774gpdeeonHHnvsW9+/+uqr11O+RkVF3XQbFApF85FSon2xDu3/PQllJRie+DWG2+9D+Pq5W1qTaTeG2xVpXYcMGUJMTAwAAwYMwGw2q7UiFQoPRBYXov3t18iP/w2Dh2P41d8Qg0e4W1azcXaV947Af4AhgAQeMJlM+5tb6X+OZHOx0Fznd6KZaV17dwrkh6NvTARVjavSulazYcMGhgwZcj3pFMCTTz6JwWBg4cKFPP7442odSYXCDcjjB/XkUFYz4u5ViGkLvP636Ozg5MvAZpPJdLvRaPQHPDtWphFaKq1rNWfOnOHFF1/kww8/vL7tlVdeITY2lrKyMh588EE+/fRT7rjjjpsTrlAonEZazEjTm/oAZM8+GH74FCK2h7tltQiNGm6j0RgOTAXuAzCZTFbgxpG6JtBQz9ib0roCZGRk8IMf/ICXX36ZXr16Xd8eG6sPdoSGhrJs2TKOHz+uDLdC0UrI9PNo//4z5GQg5iUglt3tlb7s+nCmx90HyAXeMhqNw4Bk4DGTyVTuUmUtjCvSuhYXF3Pvvffys5/9jDFjxlzfbrPZKCkpISIigqqqKrZv386UKVNatkEKheIGpGZHbklErvkAwjpiePJ5xMBb3S2rxWk0O6DRaBwNHAAmmUymg0aj8WWgxGQy/bJWuZXASgCTyTSqdvhcdnb2t/y/7mDVqlWcOnWKWbNm8atf/arBspcuXeKhhx6iqKiIoUOH8tprrxEQEMDmzZtJSUnhmWee4S9/+Qt/+9vf6NOnz/X9/vvf/xIcHMyyZcuoqqpC0zSmTJnCb37zG3x8fJzWarFY6u3ltwVUFsj2QWu22Z6bRfHLz1N18hgBE2cSvuppDGHhrVJ3NTfTXn9/f3AyO6AzhjsGOGAymXo5Pk8BfmoymRY1sJtK69oEVFrX9oFqs+vQDu9Fvvd30DTEXSsRE2a6ZQCytdK6NhoOaDKZsoArRqNxgGPTLOBUs5QpFApFCyIrK9De+Cvy9T9CbHcM//cShomzvD5qpDGcjSp5FPjAEVFyAbjfdZIUCoWiceT502hv/AXycxG3fRex6DuIJrgjWxq7JrlWbKY1HMJOGW6TyXQcGO1iLQqFQtEo0mbTU7Bu+AQiozE8/VtEv0Fu1VRktvGnLzPILrfxt4W9CPJz7dxGlWRKoVB4DTLzCtobf4X083oK1jsfQgS5dyzobF4lv9t7jVKLnadn9XO50QZluBUKhRcgNQ25Yz3y83chIADDqp8iRk10tyy2ni/iX4eziQjy5fdz4xh7S5dWGYxVhluhUHg0Mj8X7e2XIe0E3DoGw70/QnRw7xqQVrvG64ez2fZ1MSNiQ3hyUlfCA1rPv95ukky5Iq0r6FPbJ02axJQpU9i1a9f17a+//jozZsxg5syZPPzww5jNem6Wp556itmzZzN79mwefPBBysu9ah6TQtFqSCnR9u9E+/X/wsWziHt/hOFHv3C70c4tr+LZbZfZ9nUxd8RH8svp3VvVaEM7MtyuSOt69uxZ1qxZw44dO/jggw949tlnsdvtZGZm8uabb7Jx40Z27NiB3W5nzZo1ADz33HNs376d7du3061bN956660Wa6NC0VaQpSVo//w98s2/QreeGH71NwxT5ro9zO9EVjlPbrrE1WIrz07txj3Do/ExtL6mdmO4XZHWdcuWLSxdupSAgAB69uxJr169OHbsGKBPezebzdhsNiorK6+nfw0LC7teh9lsdvuNqFB4GvLEYbTnfgQphxAJ38fwkxcR0THu1SQliafy+dWOK3QI9OFPC+IY1yPMbXrc4uNOPVpBSZG9zu+am9Y1vKMPQ0bWP7rsirSuWVlZjBw58vrn6nKjR49m1apVjB07lsDAQKZNm8a0adOul3viiSfYsWMH/fv3b3TqvULRXpDmSuQnbyL3bIFucRge/zWiR293y6Kiys6rB7JIulzKxJ5hPDo+hmA/98WLQzsdnGyptK51PWCEEBQVFbFlyxYOHDhAeHg4Dz30EJ999hkrVqwA4K9//St2u51f/OIXrF27lrvvvvsmWqNQeD/y/Cm0N1+CvGw9m9/SuxF+7s/md7XEwu/2XONaiZX7RkSzbFCER7wlu8VwN9Qz9qa0rrGxsdTMyVJdbu/evfTs2ZPIyEgAFixYwJEjR64bbgAfHx+WLFnCP/7xD2W4Fe0WaatCrv0QuTlRn0zz4xcRt8S7WxYAB66U8tK+TPx8BL+e2YNbY0LcLek67abH7Yq0rnPnzuWRRx5h5cqVZGdnc/HiRUaMGIHBYODo0aNUVlYSGBjIl19+ybBhw5BScunSJXr37o2Ukm3bttGvXz+XtFeh8HTk1Uv6ZJqrFxFT5iKMDyAC3Z9Yza5JPjqRxycn8+kfGcgzU7oRHeL+3n9N2o3hjoiIYMyYMcycOZMZM2bwy1/+ssHyP//5z3n44Yf5wx/+QHx8PHfeeScAW7duJSUlhZ/85CcMGDCA2267jRkzZuDj48MLL7yAj48PI0eOZNGiRcybNw9fX1/i4+O5++67kVLy+OOPU1ZWhpSSwYMH89vf/rY1mq9QeAzSbkdu/gy57mMIDtFD/IaNdbcsAEosdv6SlMGxzHLm9O3AyjFd8PfxvBiORtO6NhOV1rUJqLSu7QPVZpDXLqO99ZI+ZX3MFH3KeivnzK6Ps3mV/GHvNQrNdh4a04W5/To2+Ritlda13fS4FQqF+5B2O3LrauTaDyAwGMOqZxCjJrlbFqAHGWw6V8QbyfrU9d/N7Un/yCB3y2oQZbgVCoVLkVlX9YiRi2dh5EQMd69ChDe9N+sKKqs0/n4wiz3pJYzqGsITE7sS1sqzIJuDMtwKhcIlSM1O+ZqP0D74F/gHIB78se4e8YBwOoArxXqoX0aplXuGRbEiPhKDh2hrDGW4FQpFiyOzM9Defpmy86dh2FgM33vE7TlGarLnUgmvHcwkwMfAczN7MMyDQv2cQRluhULRYkhNQ+7cgPz8HfD1I/yx/6MsfpTH9LKr7JK3jmaz4WwRg6KD+MnkrkQGe1aonzMow61QKFoEmZulp189exKGjsZw7yME9RtAuYdE0uSWV/GHvdc4m29m6cBO3DuiM74tnCDKRVF6N+B5AYouwlVpXRWK9o7UNLSdG/X0q1cuIu57DMOjv0R0jHS3tOsczSjjiU2XuFJs5ZkpXXlgVJcWN9rpX1vYtSULqbneeLcbw+2KtK4KRXtH5mah/fX/kB/+E/oOwvDcKxgmec4q6/osyFx+s/MqEUG+/HlBLyb2bNm4cSklZ1IrOXGkkqoqiV1r0cPXSbsx3K5I66pQtFekZkfbvgbtuUfh0jnE9x7B8PhziIhod0u7TonZxm92XeXjr/KZ3jucP86Lo1u4f4vWoWmSE0cqOXvSQo9e/sxeGIuvr+sfWk75uI1G4yWgFLADNpPJdFMrvu/Zs4fc3Nw6v2tuWtfo6GimTp1a7/euSOuqULRHZOYVtHdega/TdF/2PQ8jIqLcLetbnM6p4I9JGZSY7TwyLoY5fTu0+FuAzSY5ur+c7Awb/QcHMGBIIAaf1nnTaMrg5AyTyeQZoww3SUuldVUo2hPSZkNu+Ry5/mMICEL84EnEuGke9dvQpCTxVAHvp+TSOcSP38+Lo29EYIvXY7VoHNpbTmG+naEjg+jVP6DF62gIt0SVNNQz9qa0rgpFe0FevoD2zt/g8gXEqEmIu1Yiwj0nLht018hL+zNJzihnUs8wHhkXQ4h/y8+CrCjXOLC7jMpyjdGTgontrrtf8vLySE9PJy4ursXrrI2zhlsCW41GowT+ZTKZXq9dwGg0rgRWAphMJqKivv3qlJ2dfd310KgoJ8s1hQ4dOlBeXo6vry8dO3Zk586dDZafPHkymzZtYvny5Xz22WcsWLDAJbqqqX3sgICAG85hW8LX17dNt68uvLHNsspKuektyhPfxxDWgbCnXyRwwnSn92+tNh+/VsxzWy5QVFnFj2f0ZdnQGJe8CRTkWdi3IwO7HeYt7UZMVz2nSUZGBomJifj7+/OjH/2IgADX9sCdyg5oNBq7mkymDKPR2BnYBjxqMpn2NLCLR2YHfOSRRzh9+rRTaV3T09N5+OGHKSoqIj4+nldeecVlF0NlB2wfeFub5ddpui878wpiwkzEd36ACGnaOouubrMmJZ+fLOCDE7l0CfXj6cnd6OMC1whAXnYVh5PK8fUVjJsaSnhHvTefmZnJmjVrCAwM5IEHHkDTmhdW0pTsgE1O62o0Gp8Dykwm058aKOaRhttTUYa7feAtbZYWC3L1+8gv1kKnSAz3PIIYOqpZx3Jlm4vNNv66L5NjmeVMjtNdI65aCzLjspVjBysIDjUwflooQcF6QN7Vq1dZt24dwcHBJCQk0Lt3b89I62o0GkMAg8lkKnX8fy7wm2YpUygUHo1MO4H27quQm4WYvgCR8H1EkOd1IE5mV/CnpAxKLXb+Z2wX5vXr6LJB0gtnLZw8VklEtA9jJofg768b7fT0dNavX0+HDh1Yvnw5ISGtl+/EGadtFyDRaDRWl//QZDJtdqkqhULRqsiKMuRn7+grrEfHYPjxC4gBQ90t6wY0Kfn0ZD4fncgjJtSPX06Pc5lrRErJyeNmLp61ENPdj5HjgvFxxGhfuHCBjRs3EhERwbJly1r97bhRw20ymS4Aw262otaaw99WUOdL0RpIKeHofrSPXoeSIsScpYil9yBcPLjWHIrMNv6alMHxrAqmxoXzP+O6uMw1YrdLjh2oIPNqFb37+xM/PAjhmCJ/7tw5tmzZQnR0NEuXLiUw0DUPjoZotXBAg8FwPbxO0TA2mw2Dod1MalW4CVmQi/bhvyDlEPTsg+HRXyDiPHPx6q+yy/lzUiblVtdNqKnGatE4/GU5BXl2Bg8PpM8tAdfrSktLY9u2bcTGxnLbbbe5PHqkPlrNigYGBmI2m7FYLA2e8ICAACwWS2vJ8ghqtllKicFgcMtTXNE+kJoduWsT8vP3QNoRt9+PmL0E4eN5K7/YNYkpNQ9Taj6xYf48N6M7vTq57rdRUWbnwJ5yKss1Rk0IpmvPb6bIp6amsmPHDrp3785tt92Gn5/70sG2muEWQhAU1Pg6bt4y8t6StMc2K9yDvHoJ7b3X4MIZGDwCwz3/g4iOcbesOsktr+IvSRmcyq1keu9wVo2JIcjPdW+iRQU2Du0tR9Ng/PRQIqO/MY8pKSljR7pCAAAgAElEQVTs3r2buLg4Fi1a5HbPgfJbKBTtAFllRa43Ibd8BkEhHjldvSb7Lpfw6sEs7Bo8MTGW6b07uLS+7IwqkveV4x8gmDAjlLDwb94+kpOTSUpKok+fPsyfP9/tRhuU4VYo2jwy7QTae3+HnAx9Io3xAURoy6Y2bSksNo03knPYcr6I/pGBPDWpK7FhLZvRrzbpX1v4KrmS8I4+jJ0SQmCQ3quXUnLo0CEOHjzILbfcwpw5c/DxEHeSMtwKRRtFlpciP3kLmbRdD/F78nnEoJsOEHMZ6UUW/vTlNS4XW0kYHMFdt0bj58Jse3oebTPnTlmIjvFl9MQQfP3E9e/27dtHcnIygwYNYtasWR4VMKAMt0LRxpBSIg/vRX78bygvRSxYgVj8XYS/54X4ga5387ki3jyaQ7CfvnjviFjXTmbR7JKUIxVcvVRFz97+DB0dhMHwjdHes2cPKSkpDB06lOnTp3ucS0kZboWiDSFzMtE++hekHoVe/TE88RtEj97ullUvpRY7rxzI5ODVMkbGhvDYxFg6BrrWLFVZJUf2lZOXbWPAkED6D/4m3E9Kyc6dO0lNTWX48OFMmTLF44w2KMOtULQJZFWVnit74yfg44P47oOIGQsRBs/wydbFyewK/rwvg2KzjQdGdua2gZ0wuNhIVlZoHNpTRmmJxvCxwfTo/Y3/XNM0tm/fTlpaGqNHj2bChAkeabRBGW6FwuuRp1PQPvwnZF1DjJmiDz560EK9tbFrkv+m5vFJaj5dQv34/dxe9It0/byFkiI7B/eUYauSjJsaQnTMN3HYdrudLVu2cP78ecaPH8/YsWNdrudmUIZbofBSZEkh0vQm8uBuffDx8V8j4ke4W1aDZJdZeWlfJqdyK5nZJ5wHR7tu2npNcrOqOLKvHB8fwcSZoXTo9I3ps9lsbNq0iYsXLzJ58mRGjhzpcj03izLcCoWXITU7cs8WfeZjlUUfeFywwmMHH0H3He++VMK/DmcjZevEZldTHe4XGm5g7JRQgkO+iQ6pqqpiw4YNXL58mWnTpjFsmOdG3dREGW6FwouQ6V+jvf93uHQOBg3DcNcqREw3d8tqkDKLnX8ezmJveimDo4N4fGIsXUJdG5sN+sMi7Ssz50/r4X6jJobg5/eNz9pqtbJu3TquXbvGrFmziI+Pd7mmlkIZboXCC5CVFcg1HyB3bICwcMQPn0KMneqxg2fVHL1axK83XaSw0sY9w6JIGByJj8H1mu02ybFDFWReqSKurz9DRn4T7gdgsVhYu3YtWVlZzJs3jwEDBrhcU0uiDLdC4cFIKZFHkpD//Q+UFCKmLUAsvwcRHOpuaQ1SZZd8eCKXxFMFxIb58/t5cfSPbDxXUUtgMevZ/Qrz7QweFkifAQHfesCZzWZWr15NXl4eCxYsoF8/z8yI2BDKcCsUHorMuqbnyT51DHr2xfDIzxG9+7tbVqNcKbbw56QMLhZaWDokhrviwwn0bZ1Zh6Uldg7tKcds/vYK7NVUVFSwevVqCgsLWbRoEb17e26Me0Mow61QeBjSYkZu+C9y6xrw90d8dyVixgKPjskG/e1g49ki3j6WQ6CvgWendWPR8OavwdhUqhfz9fERTJoRSsfIb5u3srIyEhMTKS0t5bbbbqNnz56tossVKMOtUHgI190in7wJhXmIibMQK+5FhHdyt7RGKay08cqBTJIzyhnVNYRHx8fSKaj1zMvlCxZOHKkkNMzA2KnfjhwBKCkpITExkYqKCpYuXUq3bp49oNsYynArFB6AzLisu0XSTuir0az8CaLfIHfLcopDV0t59UAWlTaNlaO7sPAW1y3cW5uaiaKiuuiJovz8v113UVERiYmJWCwWli1bRmxsbKtocyXKcCsUbkRWViDXf4z8Yh0EBCLuWoWYNs/j3SIAFVV23jqaw9bzxfTuFMCTk7rSs0PrxZLb7ZLjhyrIuFxFzz7+DB317cgRgIKCAhITE7Hb7SQkJNC5c+dW0+dKnDbcRqPRBzgCXDOZTItdJ0mhaPtIKZEHdyM/fVuPFpk8B7H8e4iw1pmUcrOczK7g5QOZ5JRVtUoK1trUjBwZdGsgfQcG3NDLz8vLIzExESEEK1asIDLSc9MANJWm9LgfA04DnpmBXaHwEuTVS3oGv7Mn9Qx+jzyL6H2Lu2U5hdWu8f7xXNamFdIl1I/fzunJoM7BraqhtMTOob3lmCs1Rk0MpmuPGyfz5OTksHr1anx8fEhISKBTJ88fJ2gKThluo9HYHVgEvAA86VJFCkUbRSsvRfv438idGyA4BPG9R/Setgcl6G+Ic/mVvLQvk6slVhb078j3R3R26RqQdZGTpS8xZjAIJk4PpVPUjSYsMzOTNWvWEBgYyPLly+nQwTveYpqCsz3ul4CngbD6ChiNxpXASgCTyURUVFTzBPn6Nntfb0W1uW0j7XYqd2wg/8PXkcWFBM1bRuhdD2EI846XV5td453DV3jn0BUiQvz5y7J4xsU514Ntyet86kQRh74somOEP7MXxhIafuMq6xcvXmTNmjWEhoZy//33t7rRbq37WkgpGyxgNBoXAwtNJtPDRqNxOvBjJ3zcMiMjo1mC2uOK56rNbRd5NhXt43/DlYv4DbwV++33I+L6uluW01wutvDSvky+LjAzvZeezS80wPmB05a4zpomST1aSfrXVrp09WXk+G+WGKtJeno669evp0OHDixfvpyQENeuolMXN9Perl27Ajg1UOBMj3sSsMRoNC4EAoFwo9H4vslkuqdZ6hSKdoDMz9HXe0xOgogoxMqf0Gn+MvLz890tzSnsmmTdmQLeP55HkJ+Bn07pxoSe9b5wuwyrVSN5XwV52Tb6Dgxg0NBARB25Ti5cuMDGjRuJiIhg2bJlBAe3rt+9tWnUcJtMpp8BPwOo0eNWRluhqANpMSM3f4bckggCxG13IuYlIAJujHrwVLJKrby8X8+ZPa57KA+PjaFjK06mqaasVB+ErCjXGD42iB696w41PHfuHFu2bCE6OpqlS5cSGOj6RRncjYrjVihagOvhfZ+9A0X5eua+Fd9HRES7W5rTSCnZer6YN49mYxCCxybEMqN3uFseOLnZVSQnVSAMMGF6KJHRdZuqtLQ0tm3bRkxMDEuWLCEgwHNzkrckTTLcJpNpF7DLJUoUCi9FXjyr+7EvnIG4fhge+gmi32B3y2oS2WVWXj2YxYmsCm6NCeZ/x8cSHXLj4F9rcOm8hdSjjunrU0IIDq3bp56amsqOHTvo3r07ixcvxt/f9Tm+PQXV41Yomoksykd+/h5y/w4I74i4738RE2Z6TXgfgCYlm84W8e7xHASCh8fGMLdfB7f0sjVNcup4JRfPWekc68vICd9e+KAmKSkp7N69m7i4OBYtWoSvb/syZe2rtQpFCyCtFuT2tfqK6nYbYv4KxMI7EEHeNSCWWWrl1QOZpOZUMjw2hB+Ni3FbL7vKKkneX05ulo0+twQweFjdg5AAycnJJCUl0adPH+bPn9/ujDYow61QOI3UNOShPcjE96AgF4aPw3DHA4jO3pW0SJOSDWcKee94Lj4GwaPjY5jVxz29bIDyUjuHviynvFTj1tFBxPWt208tpeTQoUMcPHiQW265hTlz5uDj4/k5XVyBMtwKhRPIM6lon7wJ6ef1RQ0eeBwxYKi7ZTWZayVWXjmQyencSkZ1DeHhcTFEBbunlw36TMij+ypAwPjpIUR1rluLlJJ9+/aRnJzMoEGDmDVrFgYvckm1NMpwKxQNILOuon32Dhw/CJ2iEA88gRg3zav82PBNXPYHKXn4+bg3YgR0Q3zhrIVTKWbCwg2MnVz/IKSUkj179pCSksLQoUOZPn2614RWugpluBWKOpClxch1HyP3bAZff8SyexBzliL8vS/c7EqxhVcOZHImz8zY7qH8z9gYItwQl12N3SZJOVLBtfQqYrv7MXxscJ0zIUE32jt37iQ1NZXhw4czZcqUdm+0QRluheJbyCor8ot1+sCjxYyYMhex5E6vWIWmNjZNsvp0AR+fyCPQV/DkxFim9nJfLxugskJPx1pcaGfAkED6D65/YpKmaWzfvp20tDRGjx7NhAkTlNF2oAy3QoFj4PHwXn3gMT8Hbh2DYcX3EV29c13Cc/mVvHYwi4uFFib0CGPVmC5umf1Yk4I8G0eSyrHZJGMmhxDTrX7fut1uZ+vWrZw7d47x48czduzYVlTq+SjDrWj3yLMn9YHHS+egR28M338eMWiYu2U1C7NN44OUXNafKaRjoC8/m9qN8T1aP8dIbdK/tvDV0UqCgw1MmB5KWIf6o0FsNhubN2/mwoULTJ48mZEjR7aiUu9AGW5Fu0VeS0f7/F04cRg6RiLufwwxfobXDTxWczSjjH8cyianvIr5/Tty7/BoQvzdGy6n2SVfJVdw6byV6BhfRk4Ixt+//vNbVVXFhg0buHz5MtOmTWPYMO98gLoaZbgV7Q6Zn4Nc8yHywE4IDNaXDJu1BOGleS5KzDbeSM5h16USuof78+KcnsS38qo0dWExa2xZl0HWNSt9BwQw8NbAG9aErInVamXdunVcu3aNWbNmER8f34pqvQtluBXtBllWgtz4ib4CDQIxZxliwQpEqHcsaFAbKSW7L5Xwn+QcKqvsfGdoJHfER+Ln4/43huJCG4e/LMdqgRHjguneq+E8IhaLhbVr15KVlcXcuXMZOHBgKyn1TpThVrR5pMWsT1Hf8jmYzYiJMxC33YWI9J7MfbXJLrPyj0PZHMssZ0BUID8a15OeHT3jjeHKJSsnjlTg7y9YsLwbwqeswfJms5nVq1eTl5fHggUL6NevXysp9V6U4Va0WaTNhkzajlz3ERQXwrCxGJbfi+jmnZEioE+kWX+mkA9SchFCsHJ0F+b374hPAy6I1kKzS04er+TSeSuR0T6MmhhCdJdA8vLqN9wVFRWsXr2agoICFi1aRO/evVtRsfeiDLeizSGlhKP70BLfh+xr0G8QhoeeQfT3rlSrtTmTV8k/DukhfmO6hfDQGPclhaqNuVLjSFI5hfl2+twSwKBhDfuzAcrKykhMTKS0tJQlS5bQs6f3PlBbG2W4FW0GKSWcTkFb/T5cPAuxPTA88nMYNtarJ26UWey8l5LLlnNFdAry5enJXZnYM8xj2pSfayN5Xzm2KsnICcF069l4XuzS0lI+//xzKioqWLp0Kd26dWsFpW0HZbgVbQJ5/pTewz6bqq/xeN9jiAnTEQbvzR5XPfj45tEcSi12Fg/sxF23RhHs5xltklJy6ZyVk8crCQ4xMH5aKOEdG9dWVFREYmIiFouFZcuWERvrXdkVPQFluBVejUw/r/ewU4/qixl8dyVi6jyEn2e4EJrL1WIL/zyczVfZFfSPDOS5GT3oE+E5aynabJITjnwjXbr6MmJcCH7+jb8BFBQUkJiYiN1uJyEhgc6dO7eC2raHMtwKr0ReS0db8wEcOwAhYfr6jjMWe20sdjUWm8anJ/P5/FQ+Ab4GVo3pwtx+njH4WE15mZ0jSeWUFGmN5hupSV5eHomJiQghWLFiBZGRka2gtm2iDLfCq5DZGci1HyEP74HAIMSSuxCzl3jd6jN1cTSjjH8dziarrIrpvcK5f2Rnt+cXqU12RhXHDuj5s8dODaFLrHNvNjk5OaxevRofHx8SEhLo1Mn7knZ5Ep51VygU9SDzc/Q0q/t3gK8fYn4CYu5yr508U5P8iireSM4h6XIp3cL9eX5WD26NCXG3rG8hNcmZk2bOnbIQ3tGH0ZOCCaknf3ZtMjMzWbNmDYGBgSxfvpwOHTq4WG3bp1HDbTQaA4E9QICj/Kcmk+lXrhamUADIogLkRhNyz1YQIGYsQiy4HdHB+3tsVXaNtWmFmFLz0CTcPSyK5YMiPGLmY00sZo2j+yvIy7HRo5c/Q0YF4evrnOvm4sWLrF69muDgYBISEggLc3/Cq7aAMz1uCzDTZDKVGY1GP+BLo9G4yWQyHXCxNkU7RhblIzd/jtyzBTQ7YtJsxCIjIsJ7ZzvW5GhGGf8+kk1GaRXjuofywMjOxIQ1HkbX2uTn2EjeX05VlWTYmCB69nF+DCE9PZ2NGzcSFhbG8uXLCQnxrLcIb6ZRw20ymSRQPfXJz/EnXSlK0X6RhfnIzZ99Y7AnzEAsNHrdgrz1kVVq5Y2jORy6WkbXMD9+NaM7I7uGulvWDUgp+TrNQtpX5iaF+lVz4cIFNm7cSOfOnVm8eDHBwd4/BuFJCCkbt8FGo9EHSAb6Aa+ZTKZn6iizElgJYDKZRlmt1mYJ8vX1xWazNWtfb0W1Gex5OZR//h6V29eBZidwxkJCVtyLb0zbmJhhrrLz4dEM3jt8GR+D4L6xPTEO74q/r2e5RQAsZjt7v8jmyqUKevUNZdLMzg2mYq1Namoqn376KbGxsdx///34+3vem4SruJnfsuM8OeWDcspwV2M0GjsCicCjJpMptYGiMiMjw+nj1iQqKoq8vLxm7euttOc2y4Jc5KbPkF9uBSkRE2fpPuzoGHdLbBGklBy4UsYbydnkVtiYGhfOfSOjiXTjyuoNUZivz4I0myXxw4Po1c+/STM009LS2LZtGzExMSxZsoRu3bq1q3v7Zn7LXbt2BScNd5OiSkwmU5HRaNwFzAcaMtwKRYPY87LRPngd+eU23WBPmq0b7Kgu7pbWYlwttvDvI9kcz6ogrmMAry4YRI/AKnfLqpPrsyBTKgkMFEyeGUrHyKYFnaWmprJjxw66d+/O4sWL21VPu7VxJqokGqhyGO0gYDbwe5crU7RJZH4uctMn5CVtB4lusBfejohsOzPoyqx2TF/lsf5MIYG+Bh4c3ZkF/TvRpXMHj+x9Vln1WZAZV/RZkMPHBuMf0DQXTkpKCrt37yYuLo5Fixbh66sijV2JM2c3FnjH4ec2ACaTybTetbIUbQ2ZnaEPOu7fCUDQnNuwTF/s1Tmxa2PXJFvOF/HhiTzKLHZm9+3APcOj6RjouUasKN9G8v4KKis0Bt0aSN+Bzs2CrElycjJJSUn06dOH+fPnK6PdCjgTVXICGNEKWhRtEHnlInLTp8gjSeDri5g6FzF/BeG3DPLI3mdzOZpRxptHc7hSbGVIl2B+MLKzR+UWqY2UkgtnLJw+YSYwSDBxZigRUU0zuFJKDh06xMGDB+nfvz9z587Fx8czEmC1ddSjUeES5NdpaBs/0RfiDQxCzFuOmLMEEe79E2dqcrXYwptHc0jOKCcm1I+fTe3GuO6hHpNytS4sZo1jByvIzbIR292PW8cENSlqBHSjvW/fPpKTkxk0aBCzZs3C4KWLLHsjynArWgw9H/ZxtI2fwpmvIDQMsfRufbZjiOfFKt8MJRY7H3+Vx6azuh/7vhHRLB7QyeNmPdYmN6uKYwcrqLJKho4KIq5v06JGQL/Oe/fu5fjx4wwdOpTp06d79IOqLaIMt+KmkZoGxw/qPez089AxAvGdHyCmzEMEeK67oDnYNMmms4V8/FUeFVUac/t15M5bozzajw2gaZIzqWbOn7YQGt70CTXVSCnZuXMnqampDB8+nClTpiij7QY8+25TeDTSbkce3oPc+ClkXoHoGMT3HkFMmOn1+bBrI6Xk0LUy3jmWy7USK8NignlgZGd6dfL8B1NFuZ2j+ysozLfTs48/8SOczzVSE03T2L59O2lpaYwePZoJEyYoo+0mlOFWNBlprkR+uRW5bS0U5EK3OMQPn0KMnoxog4NTZ/IqeftoDqdyK+ka5s8vpnVndLcQrzBaGVespByuAHB6WbG6sNvtbN26lXPnzjF+/HjGjh3bkjIVTUQZboXTyOJC5BfrkLs3QUU59B+M4a6HYOhoRBscmMostfLe8VySLpfSIdCHVWO6MKdfR3w9aFGD+rBVSVKPVXLlopWOET6MnOB8GtYbjmWzsXnzZi5cuMDkyZMZOXJkC6tVNBVluBWNIjOvIrcmIg/sBLsdRkzAMG85os8Ad0tzCUVmG6av8th8rgg/H8F3h0aydFCEx6z12BgFeTaOHaigokKj/+AAbolvfMX1+rDZbGzYsIH09HSmTZvGsGHDWlitojkow62oEyklnDuFtjURUg6Bnz9i8hzEnKWIzl3dLc8lmG0aa08X8NmpAqx2feDxu0Oj6ORhq9DUh6ZJzp0yc/aUhaBgAxNnhBIZ3XztVquV9evXc/XqVWbNmkV8fHwLqlXcDN5xRypaDanZ4dhBtC2fw8Wzekjfbd/VQ/rC2ubKJXZN8sWFYj48kUdhpY3xPUL53vBouod7z/qVZaV2jh2ooKjATvdefgwZGYyfX/NdOhaLhbVr15KVlcXcuXMZOHBgC6pV3CzKcCsAkBYzcv8O5LY1kJOpR4jctUrP1uflC/DWR3XmvvdTcrlaYmVgVBDPTO7KoM7ekztaSsnlC1ZOHqvE4CMYNTGYrj1uLrmT2Wxm9erV5OXlsWDBAvr169dCahUthTLc7RyZn4vcuR65d6s+4Nj7Fgyr7oUR4xEG7/DpNhUpJccyy/kgJY/zBWa6h/vz06ndGO/hMx5rYzFrpBypIPuajaguenKooOCbGySuqKhg9erVFBQUsGjRInr37t1CahUtiTLc7RApJXydhty+FnlsPwBi5ETE7CXQZ4BXGa+mciqngvdTcjmZU0nnEF8emxDLtF7h+HhBpEhNsjOqSDmsz4AcPDyQPrc0PTlUbcrKykhMTKS0tJQlS5bQs2fPFlKraGmU4W5HSFsV8kgS8ot1cOkcBIcg5izT/ddtKEtfXXxdYOaDlFySM8rpFOjDytFdmNuvg8dPUa9NlVVy8lglVy5ZCevQ/BmQtSktLeXzzz+noqKCJUuW0L179xZQq3AVynC3A2RpCXLPZuTOjVBcADHdEHev0mc4trEp6bW5UmzhwxN57LtcSpi/ge8Pj2bRgE4EeOCSYY2Rk1VFyqEKLGZJ/8EB9B8ciI/Pzb8pFBUVkZiYiMViYdmyZcTGto31PdsyynC3YeS1dH3CzIFdUGWF+BEYvv8oxI9okxNmapJdZuXjr/LZdbEYfx8D3xkaydKBEYT4e5/f3lYlOZVSSfrXVkLDDYyZFNLk1Wnqo6CggMTEROx2OwkJCXTu3HYWtGjLKMPdxpA2G/LYAeSuDXD2JPj76z3rWYsRXdu+zzK7zMonqfnsuFCMQQhuG9CJFfGRdPDwJFD1kZddxfHDlVSWa/QdEMCAoS3TywbIy8sjMTERgBUrVhAZGdkix1W4Hu+8mxU3IAvzkXu2IPdugeJCiOqCuP0+fWmw0HB3y3M5WaVWPjmZz84LxQghmN+/IwnxkUR56KK8jWGzSdJOVHLxnJWQUAOTZoYScROTaWqTk5PD6tWr8fHxISEhgU6d2lae9LaOMtxejJQSzqai7dwAxw6AlDBkFIYZCyF+ZJt3h8A3BnvHhWJ8hGD+LZ1YMTjCY1dRd4a8HBsnDldQXqbRu78/A29tXja/+sjMzGTNmjUEBASQkJBAhw5tc2JVW0YZbi9EVlYgD+zUBxszr0BImD4VfdoCRHSMu+W1Cpmluktk50XdYC+8pRMJXm6wq6yS0yd0X3ZwiIEJ00OI6tKy7bl69Srr1q0jODiYhIQEwsLCWvT4itZBGW4vQl67jNy1UV9w11IJcf0Q9z2GGDMZ4d82ZzfWJrPUiilVH3RsKwYbIOtaFV8lV2A2S/rcovuyW7KXDZCens6GDRsICwsjISGBkJCQFj2+ovVo1HAbjcYewLtADKABr5tMppddLUyhIy0WZHKS7rs+fxp8/RBjpuix1737u1teq5FeZOGzk/nsTS/B1yBYdEsnlrcBg20xa6QerSTjShVhHQyMnhRCpxaKGKnJhQsX2LhxIxERESxbtozgYO+Z1q+4EWfuEBvwlMlkOmo0GsOAZKPRuM1kMp1ysbZ2jbx6Eblnqx7KV1kOXbohbr9fzx0S1vYHG6s5k1fJpyfzOXS1jEBfPUpk2eBIIrwkY199SCm5eqmKk8crsdskA4YE0m9gAIYWihipyblz59iyZQvR0dEsXbqUwMC2HbvfHmj07jeZTJlApuP/pUaj8TTQDVCGu4WRFjPy8F7kni16Zj5fP8SoiYip86B/fJueil4TKSUpWRV8ejKfr7IrCPU3cOfQKBYO6ER4gPfFYdemotzOiSOV5GbZ6BTpw7CxwYSFu6ZdaWlpbNu2jZiYGJYsWUJAG00Y1t5oUrfFaDT2AkYAB+v4biWwEsBkMhEVFdU8Qb6+zd7XW9Euncd/cyLmPVuQlRX49OhN8AOPETh9AYY22ruu6zprUrL363zePXyNtJwyIkP8eXRKb5YMiSHYCyfO1MYgfMi84sPxQ8UgYNyUKAYN7eCyB3JycjLbtm2jV69e3HXXXW4x2u3t99xa7RVSSqcKGo3GUGA38ILJZPq8keIyIyOjWYKioqLIy8tr1r7ehKwo13vXe7fqK6P7+etrNk6dB30Htvnedc3rbNMkey6V8NnJfK6WWIkJ9WNFfCQzeod7XS6R+ijItXHyuJWiAitduvoyZGQwwSGua1tKSgq7d+8mLi6ORYsW4evrHtdSe/k9V3Mz7e3atSuAUz98p66m0Wj0Az4DPnDCaCvqQWp2SDuBTPoCeeyAPg29WxxhDz5F+ZBRiOBQd0tsVcqtdraeL2LdmULyK2z06hjAU5O6MqlnmNdl66sPq0XjdIqZyxethIT6MmZyCDHdXDugmpycTFJSEn369GH+/PluM9oK1+FMVIkA3gBOm0ymv7heUttD5mQgk3YgD+yAgjw9K9+k2YhJsyCuH8HR0VS0o15JbnkVH526wJqvsqi0aQzpEszDY2MY1dU7Vk53BiklVy5aOZVixlYl6TswgAlTulFcUuDSOg8dOsTBgwfp378/c+fOxcfH+11Mihtx5lE8Cfge8JXRaDzu2PasyWTa6DpZ3o80V+gpVJO+gPOnQBj05E63P4AYPhbhd3OrlHgj5/PNrFZQ7EEAAB/nSURBVDldwJeXSxDApLhwlg2KoG9E24pyKCmy81VyBQV5djpF+XDrqGDCO/rg5+8614iUkn379pGcnMygQYOYNWsWhnYwc7a94kxUyZc46Xdp70hNg3MnkUnbkcn7wGrRU6gm3IsYPwPRqf0l8dGkJPlaOavTCkjNriDI18CSgRF8b3wffK1l7pbXolRVSc6dNHPhrAVfP8GwMUH06O3v8rcIKSV79+7l+PHjDBkyhBkzZrSZNxdF3SjnVwsgM68gD+xGHtwF+TkQFIwYNw0xaXabX1GmPiw2jV0XS1ibVsDVEiuRwb7cPzKaOX07EuLvQ1R4IHl5bcNwV8dknz5RicUs6dHbn0HDAgkIcH2PV0rJzp07SU1NZfjw4UyZMqVd3m/tDWW4m4ksKtCjQg7sgstf666QwcMQy+5BjJjQZhfYbYzsMiubzhax7esiyqwafToF8OTEWCbFhePbRgYca1KYbyP1aCVFBXY6RfowdnJQi+XKbgxN09i+fTtpaWmMHj2aCRMmKKPdTlCGuwlIcwXy6AG9Z336BEhNzxfynR/q09A7tM/UmFJKvsquYP2ZQg5f03vR43uEsXhAJwZHB7VJY2Ku1Eg7YebKJSsBgYLh44LpHufXam212+1s3bqVc+fOMX78eMaMGdMmz7OibpThbgRps8HJY8iDu5ApB8Fq1XNdL7pDd4fEtN+1+cw2jV0Xi9lwppDLxVbCAnxIGBzJ/P4diQ7x7hwi9aHZJRfPWTh70oxdg34D9SXEfP1az2jabDY2b97MhQsXmDx5MiNHjmy1uhWegTLcdSA1O5w7jTyyF3kkCcpKIDQMMXE2Yvz0duu3ria7zMpGhzuk3OEO+d/xMUzpFY5/G5kwUxspJVnXqjidYqa8TKNLV18GDw8iNKx1w+1sNhsbNmwgPT2dadOmMWzYsFatX+EZKMPtQGoaXEhDHv5SjwgpLgD/AMSwsYhx0yF+OMK3bfYincGu/f/2zjw+rqvK899aVFUqSdZqWZJlyWuseI8dmywkMQmYkBjFduAmJIQGBsL0QE9P92eYntDsDUNmYBpooOlhmrUbAhfwFhuSuN2JE4OTeIlt4sTxKtuSJWuXSqr9vdd/3KekrGiLrVJt9/v51KeqXr3lnHpVv3feufeea3GgZYAnT/Vy6OIgDgfcOKuI9y4spSFL0yFD9HTFeeVwiO5Og8JpTtbcWsCM6qn/LUSjUXbs2EFzczN33HEHixcvnnIbNOlBTgu3ZVlw9oQt1n+Ank7I88DSVWr4+bLVWT8L+nh0BmPsOtXLrlN9dIXilOW7ef+Sct69oCRjpwWbKMEBg1ePhrl4IYbX52DZ9ap7nzMFjayRSITt27fT1tbGunXraGhomHIbNOlDzgm3ZVlw/rTqEXLgD6r7ntutpvra9CE1OMaX27WKDdPipdZBnjzVy4GWASwLVlQX8PDqGayeWZg1w9FHIxoxOflKhKZTEXDANYu9zFs4tXnsRMLhMFu3bqWzs5P3vOc9zJ8/PyV2aNKHnBBuy7Kg6STWS/uUWHe0gcsFi67D0fiAEuscqxMyEt2hOP92qpenTvXSEYxT4lONjevmFzOjMPtHehqGRdOpCCdfiRCLqv7YC5f4yPenLm8fDAbZunUr3d3d3H333cyZMydltmjSh6wVbssw1CjGQ/uwDr+g0iAuFyxchuM978Ox8kYcBXq+vbhpcfDiALtP93GgZQDDgmVVfj6yspI1tUXkJaGwf7phmqquyIljYcIhi+lVbhYtz2daSWrrfAwMDLBlyxYCgQCNjY3U1dWl1B5N+pBVwm3FovDKEayX/oh15EUYCIDHA4tW4tj4kMpZF+jIGuB8b4TdZ/p45mwfvWGDEp+LxoYy1s0voWZa9kfXoO7ELl6I8dqfVE+RkjIXK97mY/okT9B7JQQCATZv3kwwGKSxsZHa2tztdqp5Mxkv3FY4iPWng3Bon3qOhNSQ82Wrcay8UeWuc7yBcYiBqMFzTf3sPtPHya4wLgdcP7OQd84rZmVNYVaObBwJy7Job41z/GiI/j6TomInq99ewIwad1r0junt7WXLli1EIhE2bNhAdXV1qk3SpBkZKdxWdyfW0f1YR/fDq4chHoeiYhxvuxXHdTdCw9Kc7rqXiGlZHG0LsvtMH89fCBA1LOqLvXx0ZSW3zZlGiS8jfwJXTOelGMf/FKany6Cg0MnKG/zU1E3diMfx6O7uZsuWLRiGwaZNm6isrEy1SZo0JCP+tZZpwrnTWEdfVCmQC2fVB9OrcKy9S4n1/AYcTl17eIhzvRH2nO3j2aZ+OoJxCjxO7phbzB3ziplf5ksboZoKLMui81KcE6+E6e4w8OWntmvfaHR2drJlyxYANm3alFNTfmneGmkr3FYkDK8exjp6QEXWfT2qkNO8Bhz3/hmOZauhelZOCdB4dAZjPNvUz56z/TT1RnA64LrqAj50XSU3zCrM2lGNo2FZFh1tcU4cUxG2L9/BkpX51M314EqzRtf29na2bt2Ky+Vi48aNlJWVpdokTRqTVsJtdXcSPPAcxh+fhuNH1dRevnwci1fC8jU4lqzCkaWT514pA1GDfecDPNPUz7FLQSzgmnIfH7++krfX514qBN7IYZ84Fqa32yDf72DpKhVhp5tgA7S2trJt2za8Xi8bN26kpKQk1SZp0py0+VdbkQjm336CQDymUiC3vhvH8jWwYJHOVw8japgcvDjInrP9HGgZIGZa1BTlcf/SCm6dPS1neoUMZ6ieyIljEfp7DfwFTpavzqe23oMzDQUboLm5mccffxy/38+mTZsoKtJdVDXjkzbC7fB6cXzkLylbeh09vkKdAhlG1DB56eIgfzgf4MXmAUJxk2Kfi3cvKGHtnGk5l7dOxDAsmpuinH4twmDApKDQyYo1+cysT68c9nDOnTvHzp07KSoqYuPGjRQW6q6qmomRNsIN4FxzK+6KChw5NHHuWIwk1kUeJzfXF/H2+mksm+HP+uHnYxGLmjSdjnL2RIRI2KK41MXKG/1U1+altWADnD17lp07d1JWVsaGDRvw+3O7zILmrZFWwq0ZX6yXzvDnTH/r0QgFTc6ciHDudAQjDtOr3Mxv8FJemR79sMfj5MmTPPnkk1RUVLBhwwZ8Pj3OQPPWGFe4hRA/AtYD7VLKJck3KfcYjMTZe66fFy4MsL/lcrG+ua6IZVUFOS/WAL3dcc6eiNByPgZATV0e8xb6KC7NnG6gx48fZ9euXVRVVdHY2Ig3R6e401wdE4m4fwJ8F/hZck3JLbpDcV5sDvDChQGOXnqNuGkxzevSYj0M07RobY5x9kSEni4Dlxtmz/cwd6EPf0FmdW88ePAgTz31FLW1taxfvx6PJzcbkTVXz7jCLaV8VggxewpsyWosy6K5P8oLzQO8cCHAia4wAFWFebx/RQ1Ly100VOTndM46kUjY5NyZKOdORQiHLAoKnSy+TnXpy0tRedWr4ciRI+zZs4f6+nruvvtu3G6dpdRcOfrXk0TipsXxjhAHWgZ4oXmAi4EoAAvKfTy4vIIbaouYVexh+vTpdOoGWUClQ5pORmk5H8U0Vf562fVeKqszI389EocOHWLv3r00NDRw++23a9HWXDWT9gsSQjwMPAwgpbzi4bputzujh/p2DkZ5vqmbfU097D/fy2DUwO10cF1tMR9YNYu3zy2jsujyvGam+3wlJPoci5qcORngtWP9dHVEcLsdXLOomGuXFlNSlrnpBMuy2LNnD3v37mXJkiXcd999qjZ8DpFrv+2p8nfShFtK+QPgB/Zb60ojyIqKioyKPg3T4kRXiIMtgxy8OMCZnggA5flubq4rZFVNIcuq/Pjz7Aa0SIDOSOCyfWSaz5NBRUUFp060ce60iq6NOBQVO1myMp/a+jzyPA7iZj+Z+rVYlsW+ffs4cOAA1157LWvXrlU1UzLVoSsk137bV+NvTU3NhNfV92xXQGcwxpHWQQ63BnmpdYBA1MTpgIaKfB5aMZ3rawqoL/Fm7K19MonFLFrORfnjv1+gqyOCywU1dR7q53koKXNlxXdmWRbPPfcchw8fZsmSJbzjHe/ICr806cNEugM+BqwFKoQQzcAXpJQ/TLZh6cRg1ODl9iBHWgc50hakuV/lqot9LlbXqqh6RVUBhd7M6ZY2lVimRUd7nOazUVpbYpgGlJZ7WLpSjW7M82SPqFmWxdNPP83LL7/MihUruOWWW7RoayadifQq+cBUGJJOxAyV/jhsC/XJrhCmBV6Xg8WVftbNL2F5lZ+6Ei9O/acclUCfwYWmKC3nooRDFnkeB7Nme5g1x8P8a2bQ1dWVahMnFdM02b17N6+++iqrVq3ipptu0qKtSQo6VYLq/XG6O8yxS0Febg9yrD1IOG7hdMD8Mh/vW1zO8qoCFlb4yMux0qhvlUjYpOV8jOamKH09Bg4HVFa7WXydhxk1ea9X58s2QTMMg6eeeoqTJ09yww03sHr16qzzUZM+5KRwxwyTk11hJdKXghzvDBGOq9b+2mkebp9bzPKqApbM8FPo0emP8YhETNqaY1y8EKOzPQ4WFJe6WHxdPjPr8vD6svtiF4/HeeKJJzhz5gw333wzq1atSrVJmiwnJ4Q7Ejd5rTPEsfYgL7eHONEZImoooa4v8XLH3GIWV/pZXOmnJD8nvpKrJhpNEOtLcSwLCgqdLLjWy8w6D0XFuXHBi8fj7Ny5k3PnznHbbbexfPnyVJukyQGyTqUsy6IzGOfVjhCvdYY43hHibE8YwwIHMKfUy7sXlLCk0s+iSj/TdIPihImETS5djNHaHKPjUhzLBH+Bk3kNXmpm5TGtJDt6hUyUaDTKjh07aG5u5vbbb2fJEl3KRzM1ZLxwxwyT090RXusMvS7W3aE4oBoTF5T72LionIaKfK6tzNepj7fIQMDgUkuMtpYY3V0GWJDvdzB3gZeaujyKS3NLrIeIRCJs376dtrY21q1bR0NDQ6pN0uQQGSXcpmXRGohxqivEqe4wr3WGOd0dJm6qtMeMwjyWzPDTUJFPw/R8Zpd4de2Pt4hlWfR2G7TZYj3QbwIwrcTJNYu8VM3Mvch6OOFwmG3bttHR0cGdd97JggULUm2SJsdIW+G2LItLAzFOdYc51RXmVLcS6WBMCYnH5WBemY/1C0tpmJ5PQ0U+pTo/fUVEwiYdbXHa22J0tMWJRiwcDiib7qZ+npeqmW78BfpOBSAYDLJ161a6u7u56667mDt3bqpN0uQgaaN0hmnxYssAF08McPRCN6e6wwxElUi7nQ7mlHq5bfY05pf7mF/mY1axjqavFNO06O0yaG+L0d4ap6/HAMDjdTC9yk1lVR6V1W483uzuDfJWGRwcZPPmzQQCARobG6mrq0u1SZocJW2E2+GAb/+xlahpUV/s4aa6IuaX5TO/3EddsZe8NJ3sNROwLItAn0lXe5zO9jid7THiMcABpeUuFi7xUVntztl89UQIBAJs3ryZYDBIY2MjtbW1qTZJk8OkjXA7HQ6+fmc9i+qrCfR2p9qcjMayLAb6TTrb43S1x+nqUOkPUL1Aamo9TK92UzHDjcejo+rx6OvrY/PmzUQiETZs2EB1dXWqTdLkOGkj3ACzir143U4C46+qScA0LPp6DXo643R3GXS1vyHU+X4HldVuKirdlFfqXPVbpaenh82bN2MYBps2baKysjLVJmk06SXcmokRCZv0dBl0d8bp6YzT22NgqjQ1+X6Vp35DqJ06/XGFdHZ2smXLFgA2bdqUU3WlNemNFu40Jxa16OtRDYi9PQa93QbBAdVo63BCSamL2fO8lFa4KKtw48vXqY/JoL29na1bt+Jyudi4cSNlZWWpNkmjeR0t3GlENGrSbwt0X7dBX4/BoC3SAD6/g5JSN7PneSgtd1Nc5nq9aJNm8mhtbWXbtm14vV42btxISUlJqk3SaC5DC3cKMAyLgX6D/j6TQJ9BJHSRzo4Q4eAb01rl+x0Ul7qpneOhpNRFcakr64s1pQMtLS1s374dv9/Ppk2bKCoqSrVJGs2b0MKdROJxi8GAwWDAJNCvRDrQp6LooakHVbrDQ3mFm6ISJdDFpS68ug/1lHP+/Hl27NhBUVERGzdupLCwMNUmaTQjooX7KjFNi+CgyWDAZMAW6YGAyWDAIBy6fGLYgkInRcUuqmflMa3YRVGxi4IiJ5WVepb3VHP27Fl27txJWVkZGzZswO/3p9okjWZUtHCPg2lahENKnEODBsFB036tnsMhi8SJu/M8DgqLnFTMcFNYpIS5sMiFv9CJ263z0enIyZMnefLJJ6moqGDDhg34fL5Um6TRjElOC7dhWERCSnzDIZNQyCSS8DoUtAgHzcuEGcCX78Bf4KR8upv8AicFhUMC7dTDxDOM48ePs2vXLqqqqmhsbMTr9abaJI1mXLJOuONxi2jEIhoxiYQTXkcsomGLcNgkEjIJhSxiUetN2ztd4PM58fkdlFW48Bfk4S9wkl/gVM/5Tpy6J0dWcOzYMXbv3k1tbS3r16/H4/Gk2iSNZkKkpXBblkU8BrGYRSxq2s/2Y4TX0YilhDliYsRH3qfTqYooeX1KhEsrnPjynfjyHfazep3ncegBKznAkSNH2LNnD/X19dx999243Wn5V9BoRmRCv1YhxJ3AtwEX8M9Sykcn2xDLsnj69wFi0X6iURPeHAxfRp7HQV6eElqP10FBkROvN88WZwcer/Oy12539k1Qq7kyDh06xN69e5k7dy533nmnFm1NxjHuL1YI4QK+B7wLaAb2CyG2SylfmUxDHA4H5dPd+AvyMYzwZcL8+rPHQV6eE3eeFmHNW8eyLPbv38/zzz/PggULWLduHS6Xrt2iyTwmEmqsAU5JKc8ACCF+CdwDTKpwA7x6ejuWZWEYxmTvOq1xuVza5ynANE16e3tpaGjgne98J06nbkjWZCYTEe6ZwIWE983A24avJIR4GHgYQEp5RQV5qqqqME0Ta3g3jizH4XBon6eIVatWccstt6REtN1ud84Vqso1n6fK34kI90g5iTf946SUPwB+MPT5lQwoWbt2LRUVFTk3GEX7PLV0d6em3rs+z9nP1fhbU1Mz4XUnEnY0A7MS3tcCF9+iTRqNRqOZJCYSce8HFggh5gAtwP3AA0m1SqPRaDSjMm7ELaWMA58CngReVYvksWQbptFoNJqRmVAHVinl74DfJdkWjUaj0UwA3R9Ko9FoMgwt3BqNRpNhaOHWaDSaDEMLt0aj0WQYjiSNXsutYYAajUYzOUyoCFOyIm7HlT6EEAevZvtMfGifc+Ohfc7+xyT4OyF0qkSj0WgyDC3cGo1Gk2Gko3D/YPxVsg7tc26gfc5+psTfZDVOajQajSZJpGPErdFoNJox0MKt0Wg0GUbKZ0m157Q8ALRIKdfb5WN/CZQBh4CHpJTRVNo42Yzg88+B64EY8CLwCSllLJU2TjbDfU5Y/h3gI1LKwpQZlyRGOM8O4CvA+wED+L6U8h9SaeNkM4LPdwBfRwWJA8CHpZSnUmnjZCKEaAICqPMZl1JeL4QoA34FzAaaACGl7JnM46ZDxP2XqHKxQ/xv4JtSygVAD/CfUmJVchnu88+BBmApkA98LBVGJZnhPiOEuB4oSY05U8Jwnz+MmpSkQUp5LSpAyTaG+/x94EEp5QrgF8BnU2JVcnmHlHKFlPJ6+/3/BHbbGrbbfj+ppFS4hRC1wN3AP9vvHcDtwG/sVX4KbEiNdclhuM+gyuZKKS0ppYWKuGtTZV8yGMlnOzL7OvA/UmVXMhnJZ+DPgS9LKU0AKWV7KmxLFqP4bAHT7NfF5MbsWfegtAuSpGGpjri/hfrjmvb7cqDXnrwB1LRpM1NhWBIZ7vPrCCHygIeAJ6baqCQzks+fArZLKVtTY1LSGcnnecB9QogDQojfCyEWpMa0pDGSzx8DfieEaEb9th9NhWFJxAKeEkIctCdMB5gx9Lu2nysn+6ApE24hxHqgXUp5MGHxSEM+s6a/4ig+J/KPwLNSyuem0KykMpLPQogaVJ73OykzLImMcZ69QNi+pf7/wI+m3LgkMYbPfwXcJaWsBX4M/P2UG5dcbpZSrgTeA3xSCHHrVBw0lRH3zUCjndz/JSpF8i2gRAgx1GiabRMTv8lnIcS/AgghvgBMB/46deYlhZHO8zFgPnDKXu4XQmRNgxWjn+dm4Lf2OluAZakxLymM5PNOYLmU8gV7nV8BN6XIvqQgpbxoP7ejzuka4JIQohrAfp70lFjKepVIKR8BHgEQQqwF/ruU8kEhxK+B96FO/p8B21Jl42Qzis8fFEJ8DHg3cMdQ/jNbGMXn9YnrCCEGpJTzU2BeUhjjPD+KunD9CLgNOJEyIyeZkXxG5XbbhBDXSClPAO9iWAN1JiOEKACcUsqA/Xod8GVgO0q7HiVJGpbqHPdI/A3w13YEVg78MMX2TAX/BMwA9gkhDgshPp9qgzRJ4VHgXiHEn4CvkZ29h17Hbqv6OPBbIcQRVI7706m1alKZAey1fXsR2CmlfAJ1nt8lhDiJulhNel5fD3nXaDSaDCMdI26NRqPRjIEWbo1Go8kwtHBrNBpNhqGFW6PRaDIMLdwajUaTYWjh1oyIEOLDQoi9V7H9T4QQX5nAenVCiAG7dknGIYR4xu6Hn4x93yKEeC0Z+9ZkNikv66qZHIQQFrAgsWSmEOKLwHwp5QdTZtg4SCnPA1Na0tUeIPKv9jDstMUufbAwmccQQnhQfcrvQ1Vq7AC2Sin/KpnH1VwdOuLWTDqZGj2nEwllH5LNI6ha8GuAIuAdwEtTdGzNFaIj7hxhKMoEvokanWoAn5FS/tj+vBxVBGgtcBx4ctj2DaiiUKtQUdnnpJTS/uwnQAioRw3lvmfYti8Dj0gpH7ff5wGtwDuBXuAskCeljAshngGeQw0NXwbsAx6QUnba234I+DtUlP4tVL32j0kp/20En+8CvoGqgd1v+/594PeAVwgxYK96DdCFqgUv7GUS+BspZcTe1z3Al4C5tv+ftEfJJR6v2v7efial/MYI9jQB/w81grAa2Ar8uZQynHB+voMqzLRLCPFDEu4MhBCzgG8Dt6CCrseklJ+yP/soalRiFWoU38NSynPDbRiB1cCWoZobqML/TRPYTpNCdMSdW1ShaiLPRAne94QQpfZn3wPCKEH5qP0AXq/JsAtVCL8S+ADwj0KIxQn7fgD4KipqG54b/xmQmK65C2iVUh4exc4HgI/Yx/Kg6l4ghFiEqqD4oG3nkC+j8UPUbEJFwBLg36WUg6hKbhellIX24yLwt8ANwApgOSoC/ax93DW2D59GpRNuZZi4CSFmA3uA744k2gk8iKpLMw91wUicWKAKNfNTPfBw4kb2XcwO4BxqZpWZ2BMxCCE2AJ8BNqEKlT0HPDaGDYk8jyox8V+EEEvtmviaNEdH3LlFDFXIP46qkTwALBRC7AfuBZbawvayEOKnKIECWA80DUXnwCEhxG9RxcCO2cu2SSn/YL8OCzEUuAIqkvycEGKalLIfFXH+yxh2/tguSoQQQgKN9vL3AY9LKffan30e+K/j+LtICHHEnjpqrOmjHgT+YmhyAyHEl1DR8edQF7kfSSl32eu2DNt2EUqAH5FSjieY35VSXrCP8VVUhD0k3ibwhYQoP3G7NUAN8OmEevVDF8hPAF+TUr5qb/e/gM8IIeonEHV/DfW9PIi6I+kSQjwipfzp2JtpUomOuLMHA8gbtiwPJV5DdCX86QGCqJTDdNRF/ELCZ4l/+HrgbUKI3qEH6o9elbBO4raXYUe0f0AVWCpBRbw/H8OXthFsBCVcrx9HShlEpThG415UdH9OCLFHCHHjGOvWcLnP5+xloFItp8fY9kGUmP9mjHWGGP4d1yS875BShkfZbhZwbtj5G6Ie+HbCuelG1bYfdxISKaUhpfyelPJm1N3EV4EfCSGunYAvmhShI+7s4TzqFjqxbOYcJlY6tAOIo8ThuL2sLuHzC8AeKeW7xtjHeNXKfoqqhucG9kkph0etE6GVhF4WQoh8VAXJEZFS7gfusXPqn0LlrWeNYutFlAAO3UHU8UYt+Auo1MZofBG4E/iFEOJ+KaUxxrqzEl4nHoNR7BriAlAnhHCPIN4XgK9KKce6GI6LlDKESp99CXUXkTUlWLMNLdzZw6+Az9olQy+iGvfeC4wVZQIq6hJCbAa+aDdyzUbVEW6yV9kBPCqEeIg3JrhdAQwM3Z5PgK2o/PQM4P9McJvh/AZ4XghxE2om8S8x8qxJQ93c3g/skFL2CSH6UXclAJeAciFEsZSyz172GOr7248S0M+jUjygcuVPCSF2AE+j8utFUsqhi1zMPtZW4F+EEB8co676J+39BFF56V9N0PcXUReuR+1JNwxglZ2e+ifg74QQh6WUx4QQxcA6KeWv7e/iGeAZKeUXR/ie/htwGHjB9uNBVDuF7lmSxuhUSfbwZeCPqLxnD0ocH5RSvjzB7T+FSkm0AT9B9TABQEoZQBWJvx91UWhD9cDwTtQ4O5r7LeouYPNEtxu2j2PAX6AuHq1AADW7SGSUTR4CmmzR/s/YDaS24D4GnLHTCzXAV1AXg6PAn4BD9jKklC+iGku/CfShGiHrh9kWRTUOVqJSDaP9t34BPAWcsR/jDlKy92+gLsTzUXdXzai+10gpt6DOxy9tX19GpaOGmIVKVY1ECPi/qHPaCXwSuFdKeWYidmlSg67HrZky7MbEayZrQJAQohDVnXCBlPLsZOwzmdjdAUfsupjEY9YCv5ZSjnvnpckcdKpEMyUIIcpQvTMeusr9vBfYjUqRfAMVHTddrX3ZipSymQmkyzSZhU6VaJKOEOLjqAa030spn73K3d2DStdcBBYA90sp9W2jJqfQqRKNRqPJMHTErdFoNBmGFm6NRqPJMLRwazQaTYahhVuj0WgyDC3cGo1Gk2H8B0LkMErHSZXyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7a6fe2001518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the call payoffs\n",
    "plt.figure(3)\n",
    "plt.plot(S_, black_scholes_call_value_six_months)\n",
    "plt.plot(S_, black_scholes_call_value_three_months)\n",
    "plt.plot(S_, black_scholes_call_value_one_month)\n",
    "plt.plot(S_, call_payoff_at_expiration)\n",
    "plt.title('Black-Schnoles price of option through time')\n",
    "plt.xlabel('Underlying stock price, S')\n",
    "plt.legend(['t=0.5', 't=0.25', 't=0.083', 't=0'], loc=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the plot above, as the time to expiration decreases, the value of the option collapses towards the value of the option at expiration, or the payoff value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Greeks\n",
    "\n",
    "A discussion about Black-Scholes would not be complete without a discussion of the Greeks. In mathematical terms, the Greeks are simply the partial derivatives of the option pricing formula. In other words, these are the sensitivities in the movement of the option price relative to the movement of other aspects of the input variables.\n",
    "\n",
    "Because the assumptions required for the Black-Scholes model do not hold true in practice, the plain vanilla formula presented above is not usually used for pricing options for the purposes of trading. (In fact if one were to rely on Black-Scholes to manage a portfolio of options, that trader would surely lose money.) It is however, largely used for understanding the sensitivities of the option with respect to it's parameters and to hedge. The most common Greeks are as follows.\n",
    "\n",
    "* **Delta**. Partial derivative of the option value with respect to the change in the underlying stock price. Delta measures how the underlying option moves with respect to moves in the underlying stock. The formula for delta is different for calls and puts.\n",
    "\n",
    "> Call delta\n",
    "> $$\\frac{\\partial V}{\\partial S} = N(d_{1})$$\n",
    "\n",
    "> Put delta\n",
    "> $$\\frac{\\partial V}{\\partial S} = N(d_{1})-1$$\n",
    "\n",
    "* **Gamma**. Second partial derivative of the option value with respect to the change in the underlying stock price. Gamma measures movements in delta or the convexity in the value of the option with respect to the underlying.\n",
    "\n",
    "> Gamma (same for call and put)\n",
    "> $$\\frac{\\partial^2 V}{\\partial S^2} = \\frac{\\phi(d_{1})}{S\\sigma \\sqrt{T-t}}$$\n",
    "\n",
    "* **Vega**. Partial derivative of the option value with respect to the change in the volatility of the underling. Vega measures how the option price moves with respect to the volatility of the underlying. (Note vega is not a real greek letter.)\n",
    "\n",
    "> Vega (same for call and put)\n",
    "> $$\\frac{\\partial V}{\\partial \\sigma} = S\\phi(d_{1}) \\sqrt{T-t}$$\n",
    "\n",
    "* **Theta**. Partial derivative of the option value with respect to the change in time. Theta measures how the value of the option decays as time passes. This was demonstrated in the plot above.\n",
    "\n",
    "> Call theta\n",
    "> $$\\frac{\\partial V}{\\partial t} = -\\frac{S\\phi(d_{1})\\sigma}{2\\sqrt{T-t}} -rKe^{-r(T-t)}N(d_2)$$\n",
    "\n",
    "> Put theta\n",
    "> $$\\frac{\\partial V}{\\partial t} = -\\frac{S\\phi(d_{1})\\sigma}{2\\sqrt{T-t}} +rKe^{-r(T-t)}N(-d_2)$$\n",
    "\n",
    "* **Rho**. Partial derivative of the option value with respect to change in the risk-free interest rate. Rho measures how the option value changes as the interest rate changes.\n",
    "\n",
    "> Call rho\n",
    "> $$\\frac{\\partial V}{\\partial r} = K(T-t)e^{-r(T-t)}N(d_2)$$\n",
    "\n",
    "> Put rho\n",
    "> $$\\frac{\\partial V}{\\partial r} = -K(T-t)e^{-r(T-t)}N(-d_2)$$\n",
    "\n",
    "Where $\\phi$\n",
    "\n",
    "> $$\\phi(x) = \\frac{e^{-\\frac{x^2}{2}}}{\\sqrt{2\\pi}}$$\n",
    "\n",
    "and $N(\\cdot)$, $d_1$ and $d_2$ are defined above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The greeks in Python\n",
    "\n",
    "Although we won't be discussing the greeks any further in this talk, let's code them up for completeness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi(x):\n",
    "    \"\"\" Phi helper function\n",
    "    \n",
    "    \"\"\"\n",
    "    return np.exp(-0.5 * x * x) / (sqrt(2.0 * pi))\n",
    "\n",
    "def call_delta(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes call delta\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: call delta\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    \n",
    "    return N(d1)\n",
    "\n",
    "def put_delta(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes put delta\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: put delta\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    \n",
    "    return N(d1) - 1.0\n",
    "\n",
    "def gamma(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes gamma\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: gamma\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    \n",
    "    return phi(d1) / (S * vol * sqrt(t))\n",
    "\n",
    "def vega(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes vega\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: vega\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    \n",
    "    return (S * phi(d1) * sqrt(t)) / 100.0\n",
    "\n",
    "def call_theta(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes call theta\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: call theta\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    d2 = d1 - (vol * np.sqrt(t))\n",
    "    \n",
    "    theta = -((S * phi(d1) * vol) / (2.0 * np.sqrt(t))) - (r * K * np.exp(-r * t) * N(d2))\n",
    "    return theta / 365.0\n",
    "\n",
    "def put_theta(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes put theta\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: put theta\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    d2 = d1 - (vol * np.sqrt(t))\n",
    "    \n",
    "    theta = -((S * phi(d1) * vol) / (2.0 * np.sqrt(t))) + (r * K * np.exp(-r * t) * N(-d2))\n",
    "    return theta / 365.0\n",
    "\n",
    "def call_rho(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes call rho\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: call rho\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    d2 = d1 - (vol * np.sqrt(t))\n",
    "    \n",
    "    rho = K * t * np.exp(-r * t) * N(d2)\n",
    "    return rho / 100.0\n",
    "\n",
    "def put_rho(S, K, r, t, vol):\n",
    "    \"\"\" Black-Scholes put rho\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :return: put rho\n",
    "    \"\"\"\n",
    "    d1 = (1.0/(vol * np.sqrt(t))) * (np.log(S/K) + (r + 0.5 * vol**2.0) * t)\n",
    "    d2 = d1 - (vol * np.sqrt(t))\n",
    "    \n",
    "    rho = -K * t * np.exp(-r * t) * N(-d2)\n",
    "    return rho / 100.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Black-Scholes call delta 0.5546\n",
      "Black-Scholes put delta -0.4454\n",
      "Black-Scholes gamma 0.0524\n",
      "Black-Scholes vega 0.1192\n",
      "Black-Scholes call theta -0.0103\n",
      "Black-Scholes put theta -0.0078\n",
      "Black-Scholes call rho 0.0978\n",
      "Black-Scholes put rho -0.1026\n"
     ]
    }
   ],
   "source": [
    "# print each of the results\n",
    "print ('Black-Scholes call delta %0.4f' % call_delta(S, K, r, t, vol))\n",
    "print ('Black-Scholes put delta %0.4f' % put_delta(S, K, r, t, vol))\n",
    "print ('Black-Scholes gamma %0.4f' % gamma(S, K, r, t, vol))\n",
    "print ('Black-Scholes vega %0.4f' % vega(S, K, r, t, vol))\n",
    "print ('Black-Scholes call theta %0.4f' % call_theta(S, K, r, t, vol))\n",
    "print ('Black-Scholes put theta %0.4f' % put_theta(S, K, r, t, vol))\n",
    "print ('Black-Scholes call rho %0.4f' % call_rho(S, K, r, t, vol))\n",
    "print ('Black-Scholes put rho %0.4f' % put_rho(S, K, r, t, vol))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implied volatility\n",
    "\n",
    "The Black-Scholes pricing formula is usually used to solve for implied volatility. Quite simply, this means setting the Black-Scholes pricing formula equal to the market observed price and using a root finding algorithm to solve for the volatility parameter which sets the difference (between model and market price) to zero.\n",
    "\n",
    "The implied volatility is generally used for making trading decisions, calibrating other more exotic securities and researching market anomolies. For example, many traders use the so called volatility skew to understand the relative value of an option to other options trading in the market. We'll see examples of the skew later.\n",
    "\n",
    "In the first of a two part [blog](http://www.bespokeoptions.com/blog/2015/10/06/implied-volatility-with-c-and-python-pt-1) series, I've discussed how to use C++ to find implied volatility. Here, we will use pure Python.\n",
    "\n",
    "First we create an objective function which we'll pass to the Scipy implementation of the [Brentq algorithm](http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.brentq.html). Then we'll create an implied volatilty function to return the volatility parameter that sets the observed market price to the model price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def call_implied_volatility_objective_function(S, K, r, t, vol, call_option_market_price):\n",
    "    \"\"\" Objective function which sets market and model prices to zero\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :param call_option_market_price: market observed option price\n",
    "    :return: error between market and model price\n",
    "    \"\"\"\n",
    "    return call_option_market_price - black_scholes_call_value(S, K, r, t, vol)\n",
    "\n",
    "def call_implied_volatility(S, K, r, t, call_option_market_price, a=-2.0, b=2.0, xtol=1e-6):\n",
    "    \"\"\" Call implied volatility function\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param call_option_market_price: market observed option price\n",
    "    :param a: lower bound for brentq method\n",
    "    :param b: upper gound for brentq method\n",
    "    :param xtol: tolerance which is considered good enough\n",
    "    :return: volatility to sets the difference between market and model price to zero\n",
    "    \n",
    "    \"\"\"\n",
    "    # avoid mirroring outer scope  \n",
    "    _S, _K, _r, _t, _call_option_market_price = S, K, r, t, call_option_market_price\n",
    "    \n",
    "    # define a nested function that takes our target param as the input\n",
    "    def fcn(vol):\n",
    "        \n",
    "        # returns the difference between market and model price at given volatility\n",
    "        return call_implied_volatility_objective_function(_S, _K, _r, _t, vol, _call_option_market_price)\n",
    "    \n",
    "    # first we try to return the results from the brentq algorithm\n",
    "    try:\n",
    "        result = brentq(fcn, a=a, b=b, xtol=xtol)\n",
    "        \n",
    "        # if the results are *too* small, sent to np.nan so we can later interpolate\n",
    "        return np.nan if result <= 1.0e-6 else result\n",
    "    \n",
    "    # if it fails then we return np.nan so we can later interpolate the results\n",
    "    except ValueError:\n",
    "        return np.nan\n",
    "\n",
    "def put_implied_volatility_objective_function(S, K, r, t, vol, put_option_market_price):\n",
    "    \"\"\" Objective function which sets market and model prices to zero\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param vol: volatility\n",
    "    :param call_option_market_price: market observed option price\n",
    "    :return: error between market and model price\n",
    "    \"\"\"\n",
    "    return put_option_market_price - black_scholes_put_value(S, K, r, t, vol)\n",
    "\n",
    "def put_implied_volatility(S, K, r, t, put_option_market_price, a=-2.0, b=2.0, xtol=1e-6):\n",
    "    \"\"\" Put implied volatility function\n",
    "    \n",
    "    :param S: underlying\n",
    "    :param K: strike price\n",
    "    :param r: rate\n",
    "    :param t: time to expiration\n",
    "    :param call_option_market_price: market observed option price\n",
    "    :param a: lower bound for brentq method\n",
    "    :param b: upper gound for brentq method\n",
    "    :param xtol: tolerance which is considered good enough\n",
    "    :return: volatility to sets the difference between market and model price to zero\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    # avoid mirroring out scope  \n",
    "    _S, _K, _r, _t, _put_option_market_price = S, K, r, t, put_option_market_price\n",
    "    \n",
    "    # define a nsted function that takes our target param as the input\n",
    "    def fcn(vol):\n",
    "        \n",
    "        # returns the difference between market and model price at given volatility\n",
    "        return put_implied_volatility_objective_function(_S, _K, _r, _t, vol, _put_option_market_price)\n",
    "    \n",
    "    # first we try to return the results from the brentq algorithm\n",
    "    try:\n",
    "        result = brentq(fcn, a=a, b=b, xtol=xtol)\n",
    "        \n",
    "        # if the results are *too* small, sent to np.nan so we can later interpolate\n",
    "        return np.nan if result <= 1.0e-6 else result\n",
    "    \n",
    "    # if it fails then we return np.nan so we can later interpolate the results\n",
    "    except ValueError:\n",
    "        return np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Call implied volatility if market and model were equal (should be close to 0.25) 0.250000\n",
      "Put implied volatility if market and model were equal (should be close to 0.25) 0.250000\n"
     ]
    }
   ],
   "source": [
    "# get the call and put values to test the implied volatility output\n",
    "call_model_price = black_scholes_call_value(S, K, r, t, vol)\n",
    "print ('Call implied volatility if market and model were equal (should be close to 0.25) %0.6f' % call_implied_volatility(S, K, r, t, call_model_price))\n",
    "\n",
    "put_model_price = black_scholes_put_value(S, K, r, t, vol)\n",
    "print ('Put implied volatility if market and model were equal (should be close to 0.25) %0.6f' % put_implied_volatility(S, K, r, t, put_model_price))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p><a name=\"fin\"></a></p>\n",
    "# Getting real options market data\n",
    "\n",
    "With the excellent [Pandas](http://pandas.pydata.org/) library (you are using Pandas, right?) we can grab entire options chains directly from Yahoo!. This is slow (it has to scrape several pages), but very useful. Chains include all options of all strikes for and expirations for a particular underlying stock.\n",
    "\n",
    "The chains are returned in a `MultiIndex pandas.DataFrame` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fix_yahoo_finance as yf\n",
    "yf.pdr_override()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a stock symbol\n",
    "#underlying_symbol = 'TSLA'\n",
    "\n",
    "# define a Options object\n",
    "# __from pandas_datareader import data as pd__\n",
    "\n",
    "#options_obj = pd.Options('TSLA', 'google')\n",
    "\n",
    "# request all chains for the underlying symbol\n",
    "# ***WARNING THIS TAKES A WHILE***\n",
    "#options_frame_live = options_obj.get_all_data()\n",
    "\n",
    "# let's pickle the dataframe so we don't have to hit the network every time\n",
    "#options_frame_live.to_pickle('options_frame.pickle')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'options_frame.pickle'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mtry_read\u001b[0;34m(path, encoding)\u001b[0m\n\u001b[1;32m     98\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 99\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mread_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mpkl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    100\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mread_wrapper\u001b[0;34m(func)\u001b[0m\n\u001b[1;32m     81\u001b[0m                             \u001b[0mcompression\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minferred_compression\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 82\u001b[0;31m                             is_text=False)\n\u001b[0m\u001b[1;32m     83\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/common.py\u001b[0m in \u001b[0;36m_get_handle\u001b[0;34m(path_or_buf, mode, encoding, compression, memory_map, is_text)\u001b[0m\n\u001b[1;32m    395\u001b[0m             \u001b[0;31m# Python 3 and binary mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 396\u001b[0;31m             \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_or_buf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    397\u001b[0m         \u001b[0mhandles\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'options_frame.pickle'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mtry_read\u001b[0;34m(path, encoding)\u001b[0m\n\u001b[1;32m    103\u001b[0m                 return read_wrapper(\n\u001b[0;32m--> 104\u001b[0;31m                     lambda f: pc.load(f, encoding=encoding, compat=False))\n\u001b[0m\u001b[1;32m    105\u001b[0m             \u001b[0;31m# compat pickle\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mread_wrapper\u001b[0;34m(func)\u001b[0m\n\u001b[1;32m     81\u001b[0m                             \u001b[0mcompression\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minferred_compression\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 82\u001b[0;31m                             is_text=False)\n\u001b[0m\u001b[1;32m     83\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/common.py\u001b[0m in \u001b[0;36m_get_handle\u001b[0;34m(path_or_buf, mode, encoding, compression, memory_map, is_text)\u001b[0m\n\u001b[1;32m    395\u001b[0m             \u001b[0;31m# Python 3 and binary mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 396\u001b[0;31m             \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_or_buf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    397\u001b[0m         \u001b[0mhandles\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'options_frame.pickle'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mread_pickle\u001b[0;34m(path, compression)\u001b[0m\n\u001b[1;32m    109\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 110\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mtry_read\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    111\u001b[0m     \u001b[0;32mexcept\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mtry_read\u001b[0;34m(path, encoding)\u001b[0m\n\u001b[1;32m    107\u001b[0m                 return read_wrapper(\n\u001b[0;32m--> 108\u001b[0;31m                     lambda f: pc.load(f, encoding=encoding, compat=True))\n\u001b[0m\u001b[1;32m    109\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mread_wrapper\u001b[0;34m(func)\u001b[0m\n\u001b[1;32m     81\u001b[0m                             \u001b[0mcompression\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minferred_compression\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 82\u001b[0;31m                             is_text=False)\n\u001b[0m\u001b[1;32m     83\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/common.py\u001b[0m in \u001b[0;36m_get_handle\u001b[0;34m(path_or_buf, mode, encoding, compression, memory_map, is_text)\u001b[0m\n\u001b[1;32m    395\u001b[0m             \u001b[0;31m# Python 3 and binary mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 396\u001b[0;31m             \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_or_buf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    397\u001b[0m         \u001b[0mhandles\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'options_frame.pickle'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mtry_read\u001b[0;34m(path, encoding)\u001b[0m\n\u001b[1;32m     98\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 99\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mread_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mpkl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    100\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mread_wrapper\u001b[0;34m(func)\u001b[0m\n\u001b[1;32m     81\u001b[0m                             \u001b[0mcompression\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minferred_compression\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 82\u001b[0;31m                             is_text=False)\n\u001b[0m\u001b[1;32m     83\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/common.py\u001b[0m in \u001b[0;36m_get_handle\u001b[0;34m(path_or_buf, mode, encoding, compression, memory_map, is_text)\u001b[0m\n\u001b[1;32m    395\u001b[0m             \u001b[0;31m# Python 3 and binary mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 396\u001b[0;31m             \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_or_buf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    397\u001b[0m         \u001b[0mhandles\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'options_frame.pickle'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mtry_read\u001b[0;34m(path, encoding)\u001b[0m\n\u001b[1;32m    103\u001b[0m                 return read_wrapper(\n\u001b[0;32m--> 104\u001b[0;31m                     lambda f: pc.load(f, encoding=encoding, compat=False))\n\u001b[0m\u001b[1;32m    105\u001b[0m             \u001b[0;31m# compat pickle\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mread_wrapper\u001b[0;34m(func)\u001b[0m\n\u001b[1;32m     81\u001b[0m                             \u001b[0mcompression\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minferred_compression\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 82\u001b[0;31m                             is_text=False)\n\u001b[0m\u001b[1;32m     83\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/common.py\u001b[0m in \u001b[0;36m_get_handle\u001b[0;34m(path_or_buf, mode, encoding, compression, memory_map, is_text)\u001b[0m\n\u001b[1;32m    395\u001b[0m             \u001b[0;31m# Python 3 and binary mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 396\u001b[0;31m             \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_or_buf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    397\u001b[0m         \u001b[0mhandles\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'options_frame.pickle'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-4a0883eccee2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# read the original frame in from cache (pickle)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0moptions_frame\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpandas\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_pickle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'options_frame.pickle'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mread_pickle\u001b[0;34m(path, compression)\u001b[0m\n\u001b[1;32m    111\u001b[0m     \u001b[0;32mexcept\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    112\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mPY3\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 113\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mtry_read\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mencoding\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'latin1'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    114\u001b[0m         \u001b[0;32mraise\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    115\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mtry_read\u001b[0;34m(path, encoding)\u001b[0m\n\u001b[1;32m    106\u001b[0m             \u001b[0;32mexcept\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    107\u001b[0m                 return read_wrapper(\n\u001b[0;32m--> 108\u001b[0;31m                     lambda f: pc.load(f, encoding=encoding, compat=True))\n\u001b[0m\u001b[1;32m    109\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    110\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mtry_read\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/pickle.py\u001b[0m in \u001b[0;36mread_wrapper\u001b[0;34m(func)\u001b[0m\n\u001b[1;32m     80\u001b[0m         f, fh = _get_handle(path, 'rb',\n\u001b[1;32m     81\u001b[0m                             \u001b[0mcompression\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minferred_compression\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 82\u001b[0;31m                             is_text=False)\n\u001b[0m\u001b[1;32m     83\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     84\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/conda3/lib/python3.6/site-packages/pandas/io/common.py\u001b[0m in \u001b[0;36m_get_handle\u001b[0;34m(path_or_buf, mode, encoding, compression, memory_map, is_text)\u001b[0m\n\u001b[1;32m    394\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    395\u001b[0m             \u001b[0;31m# Python 3 and binary mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 396\u001b[0;31m             \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_or_buf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    397\u001b[0m         \u001b[0mhandles\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    398\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'options_frame.pickle'"
     ]
    }
   ],
   "source": [
    "# read the original frame in from cache (pickle)\n",
    "options_frame = pandas.read_pickle('options_frame.pickle')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a quick look at the DataFrame that returned\n",
    "options_frame.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# and the first ten records\n",
    "options_frame.head(n=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to do some work to the `DataFrame` so that we can apply the Black-Scholes price and implied volatility to each of the options. So first let's do some clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reset the index so the strike and expiration become columns\n",
    "options_frame.reset_index(inplace=True)\n",
    "options_frame.head(n=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove PctChg and IV - we'll calculate our own\n",
    "del options_frame['PctChg']\n",
    "del options_frame['IV']\n",
    "\n",
    "# rename the columns for consistency\n",
    "columns = {'Expiry': 'Expiration',\n",
    "           'Type': 'OptionType',\n",
    "           'Symbol': 'OptionSymbol',\n",
    "           'Vol': 'Volume',\n",
    "           'Open_Int': 'OpenInterest',\n",
    "           'Underlying_Price': 'UnderlyingPrice',\n",
    "           'Quote_Time': 'QuoteDatetime',\n",
    "           'Underlying': 'UnderlyingSymbol',\n",
    "           'Chg': 'OptionChange'}\n",
    "\n",
    "options_frame.rename(columns=columns, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define some helper functions that we will apply to the `DataFrame`. These will provide some of the inputs to the option valuation functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _get_days_until_expiration(series):\n",
    "    \"\"\" Return the number of days until expiration\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: days until expiration\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    expiration = series['Expiration']\n",
    "    \n",
    "    # add the hours to the expiration date so we get the math correct\n",
    "    date_str = expiration.strftime('%Y-%m-%d') + ' 23:59:59'\n",
    "    \n",
    "    # convert date string into datetime object\n",
    "    expiry = dt.datetime.strptime(date_str, \"%Y-%m-%d %H:%M:%S\")\n",
    "    \n",
    "    # get today\n",
    "    today = dt.datetime.today()\n",
    "    \n",
    "    # return the difference and add one to count for today\n",
    "    return (expiry - today).days + 1\n",
    "\n",
    "def _get_time_fraction_until_expiration(series):\n",
    "    \"\"\" Return the fraction of a year until expiration\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: fraction of a year until expiration\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    expiration = series['Expiration']\n",
    "\n",
    "    # add the hours to the expiration date so we get the math correct\n",
    "    date_str = expiration.strftime('%Y-%m-%d') + ' 23:59:59'\n",
    "    \n",
    "    # convert date string into datetime object\n",
    "    time_tuple = time.strptime(date_str, \"%Y-%m-%d %H:%M:%S\")\n",
    "    \n",
    "    # get the number of seconds from the epoch until expiration\n",
    "    expiry_in_seconds_from_epoch = time.mktime(time_tuple)\n",
    "    \n",
    "    # get the number of seconds from the epoch to right now\n",
    "    right_now_in_seconds_from_epoch = time.time()\n",
    "    \n",
    "    # get the total number of seconds to expiration\n",
    "    seconds_until_expiration = expiry_in_seconds_from_epoch - right_now_in_seconds_from_epoch\n",
    "    \n",
    "    # seconds in year\n",
    "    seconds_in_year = 31536000.0\n",
    "    \n",
    "    # fraction of seconds to expiration to total in year, rounded\n",
    "    return max(seconds_until_expiration / seconds_in_year, 1e-10)\n",
    "\n",
    "# define terms and associated rates, these should coincide with our options\n",
    "terms = [30, 3 * 30, 6 * 30, 12 * 30, 24 * 30, 36 * 30, 60 * 30]\n",
    "rates = [0.0001, 0.0009, 0.0032, 0.0067, 0.0097, 0.0144, 0.0184]\n",
    "\n",
    "def _get_rate(series):\n",
    "    \"\"\" Interpolate rates out to 30 years\n",
    "        Note computing rates like this is not strictly theoretically \n",
    "        correct but works for illustrative purposes\n",
    "    \n",
    "    :param series: row of the dataframe, accessable by label\n",
    "    :return interpolated interest rate based on term structure\n",
    "\n",
    "    \"\"\"\n",
    "    days = series['DaysUntilExpiration']\n",
    "    \n",
    "    # generate terms for every thirty days up until our longest expiration\n",
    "    new_terms = [i for i in range(30, (60 * 30) + 1)]\n",
    "    \n",
    "    # create linear interpolation model\n",
    "    f = interp1d(terms, rates, kind='linear')\n",
    "    \n",
    "    # interpolate the values based on the new terms we created above\n",
    "    ff = f(new_terms)\n",
    "    \n",
    "    # return the interpolated rate given the days to expiration\n",
    "    return round(ff[max(days, 30) - 30], 8)\n",
    "\n",
    "def _get_mid(series):\n",
    "    \"\"\" Get the mid price between bid and ask\n",
    "    \n",
    "    :param series: row of the dataframe, accessable by label\n",
    "    :return mid price\n",
    "\n",
    "    \"\"\"\n",
    "    bid = series['Bid']\n",
    "    ask = series['Ask']\n",
    "    last = series['Last']\n",
    "    \n",
    "    # if the bid or ask doesn't exist, return 0.0\n",
    "    if np.isnan(ask) or np.isnan(bid):\n",
    "        return 0.0\n",
    "    \n",
    "    # if the bid or ask are 0.0, return the last traded price\n",
    "    elif ask == 0.0 or bid == 0.0:\n",
    "        return last\n",
    "    else:\n",
    "        return (ask + bid) / 2.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can apply these functions to each row of the `DataFrame` by setting the `axis` argument to `1` `apply` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the apply method to pass each row as a series to the various methods, returns a series in this case\n",
    "options_frame['DaysUntilExpiration'] = options_frame.apply(_get_days_until_expiration, axis=1)\n",
    "options_frame['TimeUntilExpiration'] = options_frame.apply(_get_time_fraction_until_expiration, axis=1)\n",
    "options_frame['InterestRate'] = options_frame.apply(_get_rate, axis=1)\n",
    "options_frame['Mid'] = options_frame.apply(_get_mid, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what we've done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options_frame.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "options_frame.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applying implied volatility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have all the input parameters we need to follow the same procedure to add the Black-Scholes price and the associated greeks. We'll begin as before by defining the functions we'll apply. Note that we have to compute the implied volatiltiy first which we'll use as the $vol$ input to the Black-Scholes formula. It is common practice to use the Black-Scholes formula to compute the implied volalatility then to use that result to compute the greeks.\n",
    "\n",
    "We'll take the mid price (price between bid and ask) as this is generally considered the most representative price to use. As we'll see, there will be issues with the result which we'll have to clean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _get_implied_vol_mid(series):\n",
    "    \"\"\"\n",
    "\n",
    "    \"\"\"\n",
    "    option_type = series['OptionType']\n",
    "    S = series['UnderlyingPrice']\n",
    "    K = series['Strike']\n",
    "    r = series['InterestRate']\n",
    "    t = series['TimeUntilExpiration']\n",
    "    mid = series['Mid']\n",
    "    \n",
    "    # build method name\n",
    "    meth_name = '{0}_implied_volatility'.format(option_type)\n",
    "    \n",
    "    # call from globals()\n",
    "    return float(globals().get(meth_name)(S, K, r, t, mid))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply the function to the dataframe rowwise\n",
    "options_frame['ImpliedVolatilityMid'] = options_frame.apply(_get_implied_vol_mid, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# again, let's take a look\n",
    "options_frame.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options_frame.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolating missing and bad implied volatility values\n",
    "\n",
    "In the `call_implied_volatility` and `put_implied_volatility` functions we test for the Brentq solver blowing up. The solver will blow up if the option's mid price leads to a negative implied volatility. This happens with deep in or out of the money options where market makers keep the bid-ask spread wide.\n",
    "\n",
    "Let's find the options where the implied volatility is `np.nan`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bad_iv = options_frame[np.isnan(options_frame['ImpliedVolatilityMid'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# map the count function to each strike where there is a nan implied volatility\n",
    "bad_iv.groupby(['Strike']).count()['Expiration']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll interpolate the missing implied volatility values. There are *much* more sophisticated ways of doing this (see [here](http://blog.nag.com/2013/10/implied-volatility-using-pythons-pandas.html) for example) but we'll use simple linear interpolation to prove the technique.\n",
    "\n",
    "First we'll define the function which we'll pass the frame to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _interp_implied_volatility(options_frame):\n",
    "    \"\"\" Interpolate missing (np.nan) values of implied volatility\n",
    "    We first need to split the chains into expiration and type because we cannot\n",
    "    interpolate across the entire chain, rather within these two groups\n",
    "    \n",
    "    :param options_frame: DataFrame containing options data\n",
    "    :return original DataFrame with ImpliedVolatilityMid column containing interpolated values\n",
    "\n",
    "    \"\"\"\n",
    "    # create a MultiIndex with Expiration, OptionType, the Strike as index, then sort\n",
    "    frame = options_frame.set_index(['Expiration', 'OptionType', 'Strike']).sort_index()\n",
    "    \n",
    "    # pivot the vame with ImpliedVolatilityMid as the values within the table\n",
    "    # this has Strikes along the rows and Expirations along the columns\n",
    "    # the level=1 unstack pivots on Expiration and level=0 unstack pivots on OptionType\n",
    "    unstacked = frame['ImpliedVolatilityMid'].unstack(level=1).unstack(level=0)\n",
    "    \n",
    "    # this line does three things:\n",
    "    #    first interpolates across each Expiration date down the strikes for np.nan values\n",
    "    #    second forward fills values which keeps the last interpolated value as the value to fill\n",
    "    #    third back fills values which keeps the first interpolated value as the value to fill\n",
    "    unstacked_interp = unstacked.interpolate().ffill().bfill()\n",
    "    \n",
    "    # restack into shape of original DataFrame\n",
    "    unstacked_interp_indexed = unstacked_interp.stack(level=0).stack(level=0).reset_index()\n",
    "    \n",
    "    # replace old column with the new column with interpolated and filled values\n",
    "    frame['ImpliedVolatilityMid'] = unstacked_interp_indexed.set_index(['Expiration', 'OptionType', 'Strike'])\n",
    "    \n",
    "    # give our index back\n",
    "    frame.reset_index(inplace=True)\n",
    "    \n",
    "    # return\n",
    "    return frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the completed frame\n",
    "options_frame = _interp_implied_volatility(options_frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check to see if there are any np.nans\n",
    "bad_iv_post = options_frame[np.isnan(options_frame['ImpliedVolatilityMid'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bad_iv_post.groupby(['Strike']).count()['Expiration']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applying Black-Scholes and the greeks\n",
    "\n",
    "Before we take a look at the implied volatility, let's apply the Black-Scholes formula and the associated greeks. First we define the functions we'll apply to the frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _get_option_value(series):\n",
    "    \"\"\" Return the option value given the OptionType\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: Black-Scholes option value\n",
    "\n",
    "    \"\"\"\n",
    "    option_type = series['OptionType']\n",
    "    S = series['UnderlyingPrice']\n",
    "    K = series['Strike']\n",
    "    r = series['InterestRate']\n",
    "    t = series['TimeUntilExpiration']\n",
    "    vol = series['ImpliedVolatilityMid']\n",
    "\n",
    "    meth_name = 'black_scholes_{0}_value'.format(option_type)\n",
    "    return float(globals().get(meth_name)(S, K, r, t, vol))\n",
    "\n",
    "def _get_delta(series):\n",
    "    \"\"\" Return the option delta given the OptionType\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: option delta\n",
    "\n",
    "    \"\"\"\n",
    "    option_type = series['OptionType']\n",
    "    S = series['UnderlyingPrice']\n",
    "    K = series['Strike']\n",
    "    r = series['InterestRate']\n",
    "    t = series['TimeUntilExpiration']\n",
    "    vol = series['ImpliedVolatilityMid']\n",
    "\n",
    "    meth_name = '{0}_delta'.format(option_type)\n",
    "    return float(globals().get(meth_name)(S, K, r, t, vol))\n",
    "\n",
    "def _get_gamma(series):\n",
    "    \"\"\" Return the option gamma\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: option gamma\n",
    "\n",
    "    \"\"\"\n",
    "    S = series['UnderlyingPrice']\n",
    "    K = series['Strike']\n",
    "    r = series['InterestRate']\n",
    "    t = series['TimeUntilExpiration']\n",
    "    vol = series['ImpliedVolatilityMid']\n",
    "\n",
    "    return float(gamma(S, K, r, t, vol))\n",
    "\n",
    "def _get_vega(series):\n",
    "    \"\"\" Return the option vega\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: option vega\n",
    "\n",
    "    \"\"\"\n",
    "    S = series['UnderlyingPrice']\n",
    "    K = series['Strike']\n",
    "    r = series['InterestRate']\n",
    "    t = series['TimeUntilExpiration']\n",
    "    vol = series['ImpliedVolatilityMid']\n",
    "\n",
    "    return float(vega(S, K, r, t, vol))\n",
    "\n",
    "def _get_theta(series):\n",
    "    \"\"\" Return the option theta given the OptionType\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: option theta\n",
    "\n",
    "    \"\"\"\n",
    "    option_type = series['OptionType']\n",
    "    S = series['UnderlyingPrice']\n",
    "    K = series['Strike']\n",
    "    r = series['InterestRate']\n",
    "    t = series['TimeUntilExpiration']\n",
    "    vol = series['ImpliedVolatilityMid']\n",
    "\n",
    "    meth_name = '{0}_theta'.format(option_type)\n",
    "    return float(globals().get(meth_name)(S, K, r, t, vol))\n",
    "\n",
    "def _get_rho(series):\n",
    "    \"\"\" Return the option rho given the OptionType\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: option rho\n",
    "\n",
    "    \"\"\"\n",
    "    option_type = series['OptionType']\n",
    "    S = series['UnderlyingPrice']\n",
    "    K = series['Strike']\n",
    "    r = series['InterestRate']\n",
    "    t = series['TimeUntilExpiration']\n",
    "    vol = series['ImpliedVolatilityMid']\n",
    "\n",
    "    meth_name = '{0}_rho'.format(option_type)\n",
    "    return float(globals().get(meth_name)(S, K, r, t, vol))\n",
    "\n",
    "def _get_model_error(series):\n",
    "    \"\"\" Return the error between mid price and model price\n",
    "\n",
    "    :param series: row of the dataframe, accessible by label\n",
    "    :return: error between mid price and model price\n",
    "\n",
    "    \"\"\"\n",
    "    option_mid = series['Mid']\n",
    "\n",
    "    return option_mid - _get_option_value(series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the apply method to pass each row as a series to the various methods, returns a series in this case\n",
    "options_frame['TheoreticalValue'] = options_frame.apply(_get_option_value, axis=1)\n",
    "options_frame['Delta'] = options_frame.apply(_get_delta, axis=1)\n",
    "options_frame['Gamma'] = options_frame.apply(_get_gamma, axis=1)\n",
    "options_frame['Vega'] = options_frame.apply(_get_vega, axis=1)\n",
    "options_frame['Theta'] = options_frame.apply(_get_theta, axis=1)\n",
    "options_frame['Rho'] = options_frame.apply(_get_rho, axis=1)\n",
    "options_frame['ModelError'] = options_frame.apply(_get_model_error, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options_frame.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options_frame.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyzing the model error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "It's useful to visualize some of the data we've created. A quantitative analyst has to essentially be a full stack analyst/developer understanding the implications of the valuation models and metrics. Three years ago this was called quantitative research. Now it's called data science.\n",
    "\n",
    "We'll start by simply ploting the model error. Because we computed implied volatility from the Black-Scholes formula then used that result to price the option, you may expect there to be no model error.\n",
    "\n",
    "As you can see from the histogram, this is not the case. In fact there are several occasions where the model error is relatively large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the model error\n",
    "options_frame['ModelError'].hist()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's the explanation?\n",
    "\n",
    "Based on some of my comments above, our hypothesis should be that deep in or out of the money options largely have wide bid-ask spreads leading to mis-pricing of the options. Let's plot strike v. model error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# grab the index of the 50 largest abs(errors)\n",
    "sorted_errors_idx = options_frame['ModelError'].map(abs).sort_values(ascending=False).head(50)\n",
    "\n",
    "# get the rest of the details from the frame\n",
    "errors_20_largest_by_strike = options_frame.ix[sorted_errors_idx.index]\n",
    "\n",
    "# plot model error against strike\n",
    "errors_20_largest_by_strike[['Strike', 'ModelError']].sort_values(by='Strike').plot(kind='bar', x='Strike')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also explore how the bid-ask spread contributes to model error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add a new column\n",
    "options_frame['BidAskSpread'] = options_frame['Ask'] - options_frame['Bid']\n",
    "\n",
    "# plot model error by bid-ask spread\n",
    "errors_20_largest_by_spread = options_frame.ix[sorted_errors_idx.index]\n",
    "\n",
    "# plot model error against strike, many expirations included\n",
    "errors_20_largest_by_spread[['BidAskSpread', 'ModelError']].sort_values(by='BidAskSpread').plot(kind='bar', x='BidAskSpread')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's take a look at the relationship between model error and bid-ask spread."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot a scatter plot of all errors > 1.0e-4\n",
    "options_frame[abs(options_frame['ModelError']) >= 1.0e-4].plot(kind='scatter', x='BidAskSpread', y='ModelError')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are not enough samples to get any sort of statistical relationship but we can clearly see that as the bid-ask spread increases the model error increases as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyzing implied volatility\n",
    "\n",
    "First we'll look at the so-called implied volatility skew. This plots implied volatility for one expiration across strike prices. Given all else equal, an option with a higher implied volatility will be more expensive. Generally deep in and out of the money strikes command higher prices (and by extension higher implied volatilites) than the Black-Scholes formula implies. This skew only began to appear after the crash of US stocks in October 1987. It is assumed that this is a result of assumptions in the underlying distribution assumption having \"fat tails\" rather than being completely normal.\n",
    "\n",
    "Let's take a look."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select an expiration to plot\n",
    "iv = options_frame[options_frame['Expiration'] == '2016-03-18']\n",
    "\n",
    "# get the call options\n",
    "iv_call = iv[iv['OptionType'] == 'call']\n",
    "\n",
    "# set the strike as the index so pandas plots nicely\n",
    "iv_call[['Strike', 'ImpliedVolatilityMid']].set_index('Strike').plot(title='Implied volatility skew')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The curve shows that at the money strikes have a higher implied volatiltiy out the money srikes. This is a bit unusual but very likely due to the general state of market volitlity. Traders are pricing at the money call options much more expensive than out the money call options as a result of the demand for call options at the $130 strike.\n",
    "\n",
    "In a professional setting, traders would use proprietary models (not simple linear interpolation) to rebuild the curves. If the proprietary model resulted in an implied volatility that was higher than the market, a trade would be executed to exploit the anomaly. In our case, we're simply visualizing the prevailing market state.\n",
    "\n",
    "Another interesting analysis is to look at the multi-expiration skew chart."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the monthly expirations\n",
    "expirations = options_frame['Expiration'].unique()[-5:]\n",
    "\n",
    "# get all the rows where expiration is in our list of expirations\n",
    "iv_multi = options_frame[options_frame['Expiration'].isin(expirations)]\n",
    "\n",
    "# get the call options\n",
    "iv_multi_call = iv_multi[iv_multi['OptionType'] == 'call']\n",
    "\n",
    "# pivot the data frame to put expiration dates as columns\n",
    "iv_pivoted = iv_multi_call[['DaysUntilExpiration', 'Strike', 'ImpliedVolatilityMid']].pivot(index='Strike', columns='DaysUntilExpiration', values='ImpliedVolatilityMid').dropna()\n",
    "\n",
    "# plot\n",
    "iv_pivoted.plot()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll take a look now at the term structure of volatility. This measures implied volatility at a given strike (usually an at the money strike) but across expirations. Similar to the concept of backwardization and contango in the futures markets, the term structure gives an idea of the demand for options at a particular time in the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select a strike to plot\n",
    "iv = options_frame[options_frame['Strike'] == 130.0]\n",
    "\n",
    "# get the call options\n",
    "iv_call = iv[iv['OptionType'] == 'call']\n",
    "\n",
    "# set the strike as the index so pandas plots nicely\n",
    "iv_call[['Expiration', 'ImpliedVolatilityMid']].set_index('Expiration').plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "This is a clear demonstration of a violation of the Black-Scholes model assuptions. A flat volatility term structure does not exist in the markets and we have the empirical evidence to prove it.\n",
    "\n",
    "At the $130 strike, implied is very high at very near expirations. This is a result of the recent market volatility which has increased demand for this option.\n",
    "\n",
    "One last thing to do is plot the famous volatility surface. This gives a three dimentional representation of the implied volatility against strike price and expiration. While it is sometimes hard to glean any tradeable insight from the surface directly, building surfaces is important to value other derivatives that may not have an exact strike and expiration. This usually happens in the institutional setting however so professional traders can devise their own volatility models and use the implied volatility surface to understand if options are cheap or expensive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pivot the dataframe\n",
    "iv_pivoted_surface = iv_multi_call[['DaysUntilExpiration', 'Strike', 'ImpliedVolatilityMid']].pivot(index='Strike', columns='DaysUntilExpiration', values='ImpliedVolatilityMid').dropna()\n",
    "\n",
    "# create the figure object\n",
    "fig = plt.figure()\n",
    "\n",
    "# add the subplot with projection argument\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "# get the 1d values from the pivoted dataframe\n",
    "x, y, z = iv_pivoted_surface.columns.values, iv_pivoted_surface.index.values, iv_pivoted_surface.values\n",
    "\n",
    "# return coordinate matrices from coordinate vectors\n",
    "X, Y = np.meshgrid(x, y)\n",
    "\n",
    "# set labels\n",
    "ax.set_xlabel('Days to expiration')\n",
    "ax.set_ylabel('Strike price')\n",
    "ax.set_zlabel('Implied volatility')\n",
    "ax.set_title('Implied volatility surface')\n",
    "\n",
    "# plot\n",
    "ax.plot_surface(X, Y, z, rstride=4, cstride=4, color='b')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conlusion\n",
    "\n",
    "In addition to the theory and math behind the models, we learned how to code it all in Python relying heavily on `Pandas` to gather and manipulate data and `NumPy` and `SciPy` to implement some of the models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "* http://pandas.pydata.org/pandas-docs/stable/index.html\n",
    "* http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.brentq.html\n",
    "* http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.interp1d.html\n",
    "* http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.meshgrid.html\n",
    "* http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html\n",
    "* http://nbviewer.ipython.org/github/quantopian/qgrid/blob/master/qgrid_demo.ipynb\n",
    "* http://daringfireball.net/projects/markdown/syntax\n",
    "* http://www.codecogs.com/latex/eqneditor.php\n",
    "* http://www.bespokeoptions.com/blog/2015/10/06/implied-volatility-with-c-and-python-pt-1\n",
    "* http://blog.nag.com/2013/10/implied-volatility-using-pythons-pandas.html\n",
    "* http://www.888options.com\n",
    "* https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model\n",
    "* https://en.wikipedia.org/wiki/Greeks_(finance)\n",
    "* https://en.wikipedia.org/wiki/Volatility_smile\n",
    "* http://pythonpodcast.com/scott-sanderson-algorithmic-trading.html\n",
    "* http://pythonpodcast.com/yves-hilpisch-quant-finance.html"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}