4.1 Option pricing with transaction costs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Davis-Panas-Zariphopoulou model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook presents a model for pricing options in a market with proportional transaction costs. The model is taken from the celebrated work of Davis-Panas-Zariphopoulou 1993 [*link-to-the-paper*](https://web.ma.utexas.edu/users/zariphop/pdfs/TZ-7.pdf). \n",
    "\n",
    "#### This is a very powerful model!\n",
    "\n",
    "However, due to its complexity (and time complexity), it is not very well known to the practitioners.\n",
    "\n",
    "The purpose of this notebook is to explain **in simple terms** the main ideas of the model, and show how to implement it numerically. **The results will surprise you!**\n",
    "\n",
    "## Contents\n",
    "   - [Model description](#sec1)\n",
    "      - [Portfolio dynamics (original)](#sec1.1)\n",
    "      - [Some definitions](#sec1.2)\n",
    "   - [Singular control problem](#sec2)\n",
    "      - [Maximization problem](#sec2.1)\n",
    "      - [Indifference pricing](#sec2.2)\n",
    "   - [Variable reduction](#sec3)\n",
    "      - [Minimization problem](#sec3.1)\n",
    "      - [Portfolio dynamics (2 state variables)](#sec3.2)\n",
    "      - [HJB variational inequality](#sec3.3)\n",
    "      - [Indifference price (explicit form)](#sec3.4)\n",
    "   - [Numerical Solution](#sec4)\n",
    "       - [Discrete SDE](#sec4.1)\n",
    "       - [Algorithm](#sec4.2)\n",
    "   - [Numerical Computations](#sec5)\n",
    "       - [Time complexity](#sec5.1)\n",
    "       - [Is the risk aversion important?](#sec5.2)\n",
    "       - [Is the drift important?](#sec5.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from IPython.display import display\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec1'></a>\n",
    "# Model Description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec1.1'></a>\n",
    "### Portfolio dynamics (original)\n",
    "Let us consider a portfolio composed by: \n",
    "- A bank account $\\mathbf{B}$ paying an interest rate $r > 0$. \n",
    "- A stock $\\mathbf{S}$. \n",
    "- A number of shares $\\mathbf{Y}$ of the stock $\\mathbf{S}$. \n",
    "\n",
    "The 3-D state of the portfolio at time $t\\in [t_0,T]$ is $(B_t,Y_t,S_t)$ and evolves following the SDE:\n",
    "\n",
    "\\begin{equation}\n",
    " \\begin{cases}\n",
    " dB_t &=  rB_t dt - (1+\\theta_b)S_{t} dL_t + (1-\\theta_s) S_{t} dM_t \\\\\n",
    " dY_t &=  dL_t - dM_t \\\\\n",
    " dS_t &=  S_{t} \\left( \\mu dt + \\sigma dW_t \\right).\n",
    "\\end{cases}\n",
    "\\end{equation} \n",
    "\n",
    "The parameters $\\theta_b$, $\\theta_s \\geq 0$ are the proportional transaction costs when buying and selling respectively. \n",
    "\n",
    "The processes  $\\{(L_t, M_t)\\}_{t \\in [t_0,T]}$ are the **trading strategies**, i.e. the **controls** of the problem. \n",
    "\n",
    "The process $\\{L_t\\}_{t}$ represents the cumulative number of shares bought up to time $t$, and $\\{M_t\\}_{t}$ represents the number of shares sold up to time $t$. \n",
    "They are right-continuous, finite variation, non-negative and increasing processes. If the time $t$ is a discontinuous point (there is a transaction), the variation of the processes are indicated as\n",
    "$$ \\Delta L_t= L(t)-L(t^-) \\quad \\quad \\Delta M_t= M(t)-M(t^-) $$\n",
    "\n",
    "Let us consider an example. If at time $t$, the investor wants to buy $\\Delta L_t$ shares. Then the portfolio changes as\n",
    "\n",
    "$$ \\Delta Y_t =  \\underbrace{\\Delta L_t}_{\\text{shares bought}} \\quad \\quad\n",
    " \\Delta B_t =  - \\underbrace{(1+\\theta_b)S_t}_{\\text{adjusted price}} \\Delta L_t $$\n",
    "\n",
    "where the **adjusted price** is the real cost of the stock (incorporating the transaction cost).\n",
    "\n",
    "If there are no transactions, the portfolio has the simple well known evolution:\n",
    "\n",
    "$$\n",
    " \\begin{cases}\n",
    " dB_t &=  rB_t dt \\\\\n",
    " dY_t &=  0 \\\\\n",
    " dS_t &=  S_{t} \\left( \\mu dt + \\sigma dW_t \\right).\n",
    "\\end{cases}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec1.2'></a>\n",
    "### Some definitions\n",
    "\n",
    "The **cash value** function $c(y,s) : \\mathbb{R} \\times \\mathbb{R}^+ \\to \\mathbb{R}$, is defined as the value in cash when the shares in the portfolio are liquidated i.e.  \n",
    "long positions are sold and short positions are covered.\n",
    "\n",
    "$$\n",
    "c(y,s) := \\begin{cases} \n",
    "(1+\\theta_b)ys, & \\mbox{if } y\\leq 0 \\\\ \n",
    "(1-\\theta_s)ys, & \\mbox{if } y>0 . \n",
    "\\end{cases} \n",
    "$$\n",
    "\n",
    "For $t\\in [t_0,T]$, the **total wealth** process in a portfolio with zero options is defined as:\n",
    "\n",
    "$$ \\mathcal{W}^0_t := B_t + c(Y_t,S_t). $$\n",
    "\n",
    "If the portfolio contains an option with maturity $T$ and strike $K$, then the wealth process becomes:    \n",
    "**Writer**:\n",
    " \n",
    "$$ \\mathcal{W}^{w}_t = \\; B_t + c(Y_t,S_t) \\mathbb{1}_{\\{t < T\\}} +\n",
    " c(Y_t,S_t) \\mathbb{1}_{\\{t = T,\\, S_t(1+ \\theta_b ) \\leq K\\}} +\n",
    " \\biggl( c\\bigl( Y_t-1,S_t \\bigr) + K \\biggr) \\mathbb{1}_{\\{t=T,\\, S_t(1+ \\theta_b ) > K \\}} $$\n",
    " \n",
    "**Buyer**:   \n",
    "\n",
    " $$   \\mathcal{W}^{b}_t = \\; B_t + c(Y_t,S_t) \\mathbb{1}_{\\{t < T\\}} +\n",
    "  c(Y_t,S_t) \\mathbb{1}_{\\{t = T,\\, S_t(1+ \\theta_b ) \\leq K\\}} +\n",
    "  \\biggl( c\\bigl( Y_t+1,S_t \\bigr) - K \\biggr) \\mathbb{1}_{\\{t=T,\\, S_t(1+ \\theta_b ) > K \\}}  $$\n",
    "  \n",
    "For $t_0 \\leq t<T$, the wealths $\\mathcal{W}^{w}_t$ and $\\mathcal{W}^{b}_t$ are equal to $\\mathcal{W}^{0}_t$, but when $t = T$ they differ because of the payoff of the option. The writer gives away a share and recives the strike and the buyer receive a share and pays the strike. \n",
    "\n",
    "Note that considering a market with transaction costs, implies a different condition for the exercise of the option. Now the buyer should exercise if $S_t(1+ \\theta_b ) > K$, because the true price of the share incorporates the value of the transaction costs. Let's see the plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcjXX7wPHPZd+z9iRLVAhjjDFJm6KSIlotUfKztC/apEXrozyK0qZFSoTKkiyJQupJwoOyZIsMlX3fZrh+f3zvmY4xY8b4njmzXO/X67ycc5/7vs733HOc69zfVVQVY4wxBiBfpAtgjDEm+7CkYIwxJpklBWOMMcksKRhjjElmScEYY0wySwrGGGOSWVIweYKI1BKR/4nIbhG5T0SKisiXIrJTRD6LdPl8EZElInJppMthcq4CkS6AyXtEZC3wL+AwsBeYDNyrqnvC+LKPAjNVtUFQhluCMpRT1cQwvm6WUtW6kS6DydnsSsFEyjWqWgKIBc4Fngzz650BLEnxeEVuSQgiYj/wjBeWFExEqeoGYAoQBSAiXURkWVDNs0ZEbk/aV0R+FZFrQh4XFJEtIhITPG4dVJ/sEJGZIlI72P4t0BR4Q0T2iMhIoA/QLnjcNWW5RCS/iDwuIquDsswXkSrBcxeIyM9B1dPPInJByHEzReQFEflvEPtLESknIiNEZFewf7WQ/TWozloTvJf+IpIveO4sEflWRLYGz40QkdIhx64VkV4ishjYKyIFgm2XB883EpF5wev+LSIDQo5N9VyFxH1YRBYH73G0iBTJ1B/Y5Dyqaje7ZekNWAtcHtyvgvsF/3zwuCVwFiDAJcA+IDZ47lFgdEicNsAvwf2auKqoK4CCwb6rgELB8zOBbiHHPgMMP04ZHwF+AWoFZakPlAPKAtuBW3DVrx2Cx+VCXmdV8B5OAZYCK4DLg/2HAUNDXkeBGUHcqsG+3YLnzg7eT2GgAvAd8GqK87gwOIdFUzm3PwK3BPdLAI0zeK7WAnOB04NyLQPuiPTnxm5Zc7MrBRMp40VkB/A9MAvoC6Cqk1R1tTqzgK+Bi4NjhgNXi0ip4PEtwMfB/XbAJFWdpqoJwMtAUSD5V/wJ6gY8qaq/BWVZpKpbcUlrpap+rKqJqjoSWA5cE3Ls0OA97MRdBa1W1enqqqo+AxqkeK1+qrpNVf8AXsUlGlR1VfB+DqrqZmAALlGGGqSq61V1fyrvIQE4W0TKq+oeVZ0TbM/IuRqkqhtVdRvwJRCT8VNncjJLCiZSrlXV0qp6hqrelfSlJiJXicgcEdkWJI2rgfIAqroR+AG4IahGuQoYEcQ7HViXFFxVjwDrgUqZLF8VYHUq2496ncC6FK/zd8j9/ak8LpHi+PUpYp0OICKnisgoEdkgIrtwSbH8cY5NqSvuqmB5UG3VKrX3kMa5+ivk/r5UymxyKUsKJtsQkcLAGNwv13+pamlczyQJ2e0joBNwE/CjujYJgI24xuOkWIL7Yt9A5qzHVQGldNTrBKqexOuAK2dorI3B/Rdx1UvRqloK974lxbFpTnOsqitVtQNwKtAP+FxEiuP/XJlcxJKCyU4K4erPNwOJInIV0DzFPuNxPZbux9XPJ/kUaCkil4lIQeAh4CDw30yW5X3geRGpIU60iJTDJamaInJz0LDbDqgDTMzk6wA8IiJlgobs+4HRwfaSwB5gh4hUwrVzZJiIdBKRCsGVwI5g82H8nyuTi1hSMNmGqu4G7sN9aW0HbgYmpNhnP+5qojowNmT7b7hf0q8DW3B1/Neo6qFMFmdAUI6vgV3AEFxj7lagFe6LdCuukbaVqm7J5OsAfAHMxzUaTwpeC+BZXALcGWwfm+rRaWsBLBGRPcBrQHtVPRCGc2VyEVG1RXZMziIifYCaqtop0mU5WSKiQA1VXRXpshgDNqLZ5DAiUhbXgHpLpMtiTG5k1UcmxxCR7rgG4Cmq+l2ky2NMbmTVR8YYY5LZlYIxxphkYWtTEJEPcL00NqlqVCrPdwR6BQ/3AHeq6qL04pYvX16rVavms6jGGJPrzZ8/f4uqVkhvv3A2NH8IvMHRfclD/Q5coqrbg/7o7wLnpRe0WrVqzJs3z1shjTEmLxCRlCPxUxW2pKCq34XOBpnK86EDZeYAlcNVFmOMMRmTXdoUuuImDkuViPQIpgCet3nz5iwsljHG5C0RTwoi0hSXFHqltY+qvquqcaoaV6FCulVixhhjMimig9dEJBo3x8xVwfQBmZKQkEB8fDwHDhzwVzhjjqNIkSJUrlyZggULRrooxngVsaQgIlVxc7ncoqorTiZWfHw8JUuWpFq1argJH40JH1Vl69atxMfHU7169UgXxxivwtkldSRwKVBeROKBp3GrPKGqg3HLIZYD3gq+yBNVNS4zr3XgwAFLCCbLiAjlypXD2rdMbhTO3kcd0nm+G251Ky8sIZisZJ83k1tFvKHZGGNMOlTh2WdhUbrje0+aJQUPxo0bR0xMzFG3fPnyMWVKmr1sI2bHjh289dZbaT6fP39+YmJiiIqK4qabbmLfvn1eX3/QoEHUrl2bjh07cvDgQS6//HJiYmIYPXp0+gcbk1cNGADPPAOjRoX/tVQ1R90aNmyoKS1duvSYbZH0zjvvaJMmTfTw4cMZ2v/IkSMZ3vdk/f7771q3bt00ny9evHjy/ZtvvllfeeUVr69fq1YtXbNmjaqq/vjjj9qkSROv8bNSdvvcmVxq5EhVUL3pJtWT+J4A5mkGvmPtSsGzFStW8Nxzz/Hxxx+TL587vf379+fcc88lOjqap59+GoC1a9dSu3Zt7rrrLmJjY1m/fj0jR46kXr16REVF0atX6sM2fv75Zy644ALq169Po0aN2L17NwcOHKBLly7Uq1ePBg0aMGPGDACWLFlCo0aNiImJITo6mpUrV/LYY4+xevVqYmJieOSR46/uePHFF7NqlVv75dprr6Vhw4bUrVuXd999F4AhQ4bQs2fP5P3fe+89HnzwQQAGDBhAVFQUUVFRvPrqqwDccccdrFmzhtatW9OvXz86derEwoULiYmJYfXq1Zk95cbkXjNnQufOcPHFMGwY5MuCr+yMZI7sdEvvSuH++1UvucTv7f77M5CGVfXQoUPasGFDHTlyZPK2qVOnavfu3ZOvBlq2bKmzZs3S33//XUVEf/zxR1VV3bBhg1apUkU3bdqkCQkJ2rRpUx03btxR8Q8ePKjVq1fXuXPnqqrqzp07NSEhQV9++WW97bbbVFV12bJlWqVKFd2/f7/ec889Onz48ORj9+3bl+ErhYSEBG3durW+9dZbqqq6detWVVXdt2+f1q1bV7ds2aJ79uzRM888Uw8dOqSqqueff74uXrxY582bp1FRUbpnzx7dvXu31qlTRxcsWKCqqmeccYZu3rxZVVVnzJihLVu2zNjJzYbsSsGE1a+/qp5yimrt2qrB/7+TgV0pZL2nnnqKunXr0r59++RtX3/9NV9//TUNGjQgNjaW5cuXs3LlSgDOOOMMGjduDLgrgEsvvZQKFSpQoEABOnbsyHffHb2OzG+//UbFihU599xzAShVqhQFChTg+++/55Zb3EJk55xzDmeccQYrVqzg/PPPp2/fvvTr149169ZRtGjRdN/D/v37iYmJIS4ujqpVq9K1a1fAtQXUr1+fxo0bs379elauXEnx4sVp1qwZEydOZPny5SQkJFCvXj2+//57rrvuOooXL06JEiW4/vrrmT179smfYGPyig0b4KqroFgxmDIFypbNspfOdctxBjUVWW7mzJmMGTOGBQsWHLVdVenduze33377UdvXrl1L8eLFj9ovPaqaalfItI69+eabOe+885g0aRJXXnkl77//PmeeeeZxX6No0aIsXLjwqG0zZ85k+vTp/PjjjxQrVoxLL700efR4t27d6Nu3L+eccw5dunTJ8HsxxqRh5064+mrYvh1mz4YzzsjSl7crBQ+2b99Oly5dGDZsGCVLljzquSuvvJIPPviAPXv2ALBhwwY2bdp0TIzzzjuPWbNmsWXLFg4fPszIkSO55JJLjtrnnHPOYePGjfz8888A7N69m8TERJo0acKIESMA16bxxx9/UKtWLdasWcOZZ57JfffdR+vWrVm8eDElS5Zk9+7dJ/T+du7cSZkyZShWrBjLly9nzpw5R5V7/fr1fPLJJ3To4IamNGnShPHjx7Nv3z727t3LuHHjuPjii0/oNY3Jkw4dghtugKVLYexYiInJ8iLkuiuFSBg8eDCbNm3izjvvPGp77969adeuHcuWLeP8888HoESJEgwfPpz8+fMftW/FihV58cUXadq0KarK1VdfTZs2bY7ap1ChQowePZp7772X/fv3U7RoUaZPn85dd93FHXfcQb169ShQoAAffvghhQsXZvTo0QwfPpyCBQty2mmn0adPH8qWLcuFF15IVFQUV111Ff3790/3/bVo0YLBgwcTHR1NrVq1kqu8krRt25aFCxdSpkwZAGJjY7ntttto1KgR4K4mGjRocGIn1Zi8RhW6doVvvoEPP4QrrohIMXLcGs1xcXGacpGdZcuWUbt27QiVyLRq1YqePXty2WWXRbooWco+d8arxx+HF1+EF16AJ57wHl5E5msGphKy6iOTaTt27KBmzZoULVo0zyUEY7x6+22XEHr0cMkhgqz6yGRa6dKlWbHipCa4NcZMmAD33AOtWsGbb0KE59WyKwVjjImUn36C9u2hYUM3hUWByP9Ot6RgjDGRsHKluzo4/XSYOBFCuqhHkiUFY4zJaps2ucFpAF99BaeeGtnyhIj8tYoxxuQle/e6K4SNG2HGDDj77EiX6Ch2peDB1q1bk6fMPu2006hUqVLy40OHDkW6eIwdO5bly5cnP37iiSeSJ82LtJRlS8/06dM55ZRTiImJoXbt2vz73//2Wh5VpW3btkRHRzNo0CCWLl1K/fr1adCgAWvXrvX6WiYPSkx0bQjz57s2hPPOi3SJjmFXCh6UK1cueWqIZ555hhIlSvDwww8ftU/yZFNZMcthCmPHjiVfvnycc845AN6/SE9GyrJlRNOmTRk/fjx79uwhOjqaVq1aUb9+fS/l2bBhA/Pnz0+etfWFF17gxhtv5KmnnvIS3+RhqnD33a794O23oXXrSJcoVXalEEarVq0iKiqKO+64g9jYWP7880969OhBXFwcdevW5bnnnkvet3LlyjzzzDM0aNCA6Ojo5K6e3377LfXr1ycmJobY2Fj27t3Lrl27aNasGbGxsURHRzNx4sTkOEOHDiU6Opr69evTpUsXZs+ezeTJk+nZsycxMTGsXbuWTp06MX78eACmTZtGTEwM9erVo3v37slXNmmVJ1RiYiI9e/YkKiqK6Ojo5MV70or5yCOPUKdOHaKjo+nVq1eqZRs4cCB16tShfv36dOrU6bjnt0SJEsTGxrJ69WpWr17NxRdfTIMGDWjYsCE//fQTAB06dGDSpEnJx7Rr147Jkyezf/9+OnfuTL169YiNjU2efLB58+Zs3LiRmJgYnnvuOd544w0GDx7M5ZdffmJ/fGNS6tsX3n0XeveGO+6IdGnSlpGpVLPTLd1FdiI5d7aqPv3009q/f39VVV25cqWKSPJU16r/TEGdkJCgF110kS5ZskRVVStVqpQ8TfVrr72mt99+u6qqtmjRQufMmaOqqrt379bExEQ9dOiQ7tq1S1VV//77bz377LNVVXXhwoVaq1at5NdI+rdjx45HTcOd9Hjv3r1auXJlXbVqlaq6RXVef/3145Yn1KBBg7Rt27aamJiY/Hppxfzrr7+0Tp06euTIEVVV3b59e6plO+200/TgwYNH7RNq2rRp2qZNG1VV3bRpk1apUkWXL1+ue/fu1f3796uqmz68UaNGqqo6ffp0veGGG1RVddu2bVq9enVNTEzUl156Sbt166aqqr/++qtWrVpVDx48qCtXrtT69esnv94TTzyhAwcOPKYcqjZ1tjkBH37oFsrp1Ek1+D+Q1bCps7OHs846K3mqa4CRI0cSGxtLbGwsy5YtY+nSpcnPXX/99QA0bNgwuf76wgsv5IEHHuD1119n165d5M+fH1WlV69eREdH07x5c9avX8+WLVv49ttvadeuHWWDaXbLpjPd7rJly6hRowZnnXUWALfeeutR03WnVp5Q06dP54477kiex6ls2bJpxixbtiz58uWje/fujBs37qgZYkPVrVuXTp06MWLECAoWLJjqPjNmzKBBgwa0aNGCp556ilq1anHw4EG6du1KVFQU7du3Tz6vzZo1Y+nSpWzdupURI0bQtm1b8ufPf9R043Xr1uX0009PXlDIGK++/hq6dYPLLoMhQyI+OC09ua9NIVJzZ6ch9Mtv5cqVvPbaa8ydO5fSpUvTqVOn5CmoAQoXLgy4dZITExMBePLJJ2ndujWTJk3i3HPPZebMmcyaNYudO3eyYMECChQoQOXKlTlw4ECaU2unRdOZ9yq18qQ8PuXrpRWzYMGCzJs3j2nTpjFq1Cjefvttvv7662P2mzp1KrNmzeKLL77ghRde4Ndffz1m8sCkNoVQr7zyClWqVGH48OEkJCRQokQJAESEjh078sknn/Dhhx/yySefZOi9G+PFwoVu1tM6dWDMGChUKNIlSpddKWShXbt2UbJkSUqVKsWff/7J1KlT0z1m9erVREdH07t3bxo0aMBvv/3Gzp07OfXUUylQoADTpk1jw4YNAFx++eWMGjWKbdu2AST/m9Z02XXq1GHlypWsWbMGgOHDhx8zXffxNG/enLfffpvDhw8nv15aMXfv3s2uXbto1aoVAwcO5H//+98xZTt8+DDx8fE0a9aM/v37s3nzZvbt25ehsuzcuZOKFSsiInz00UdHfel36dKF/v37U6RIEWrVqgVw1HTjy5Yt488//+TsbNY10ORw69a5dRFKl4bJk+GUUyJdogyxpJCFYmNjqVOnDlFRUXTv3p0LL7ww3WNefvnl5Ibc0qVL07x5c2655Rb++9//EhcXx2effUaNGjUAiI6O5tFHH6VJkyZHrcHcoUMH+vbtm9yYm6RYsWIMGTKE66+/nnr16lG4cGG6d++e4fdz++23c9pppyU3bH/66adpxty5cyctW7akfv36NGvWjAEDBhxTtlWrVnHzzTcTHR1NbGwsvXr1OmZ9irTcc889vP/++zRu3Jh169YlX+UAnH766dSsWTN5ESAgefrxevXq0bFjR4YNG0ahHPArzuQQ27a5wWn79rnBaZUqRbpEGWZTZ5tcb+/evdSrV49FixZlOMlkhH3uTKoOHIDmzd28RlOnwqWXRrpEQDaYOltEPhCRTSLyaxrPi4gMEpFVIrJYRGLDVRaTd02dOpXatWvTs2dPrwnBmFQdOQKdO7tlND/6KNskhBMRzobmD4E3gGFpPH8VUCO4nQe8HfxrjDdXXnklf/zxR6SLYfKKRx6BTz+F/v3dyOUcKGxXCqr6HbDtOLu0AYYFXWjnAKVFpOJJvF5mDzXmhNnnzRzj1VdhwAC47z546CHv4b/9Fnbt8h72GJFsaK4ErA95HB9sO2FFihRh69at9h/VZAlVZevWrRQpUiTSRTHZxeefw4MPuu6nAwZ4H4swc6Zrt370Ua9hUxXJcQqpnbVUv9VFpAfQA6Bq1arHPF+5cmXi4+PZvHmz1wIak5YiRYpQuXLlSBfDZAezZ0OnTnDBBfDxx5BiXM3J+uUXuPZaOOssN1NGuEUyKcQDVUIeVwY2prajqr4LvAuu91HK5wsWLEj16tXDUUZjjEnbsmXQpg1Ur+6W1Sxa1Gv49evdFULx4q5nazqTFHgRyeqjCcCtQS+kxsBOVf0zguUxxpiM27gRWrSAwoVhyhTv39jbtrnwu3e7hJBKJUlYhO1KQURGApcC5UUkHngaKAigqoOBycDVwCpgH9Al9UjGGJPN7NoFLVvC1q0waxZUq+Y1/P797gJk1SqXEOrV8xr+uMKWFFS1QzrPK3B3uF7fGGPCIiEBbrzRVfZPnAgNG3oNf/gwdOwI338Po0dD06Zew6cr902IZ4wx4aLqZjydNg0++MDV73gOf//9MG4cDBwIbdt6DZ8hNveRMcZkVJ8+MGwYPPssdPFf4/3SS/Dmm/Dww/DAA97DZ4glBWOMyYh33oEXXoCuXSEMy7N+9BE8/jh06AD9+nkPn2GWFIwxJj0TJ8Jdd7mpsAcP9j447auv/lmH58MPIQJLuSezpGCMMcczdy60awcNGriW3wJ+m2Lnz3ft1lFRMHZs5NfhsaRgjDFpWbUKWrWCf/0LJk2CYEU/X1avdhcf5cu7dXhKlfIaPlMsKRhjTGo2b3bDiY8ccfU7//qX9/AtWkBiolt2oWKmpwP1y7qkGmNMSvv2wTXXQHy8m560Zk2v4ffudWPfNmyAb76BYJXYbMGSgjHGhEpMdGshzJ3rKvnPP99r+IQEN/5g/nw3HsFz+JNmScEYY5KouvUQvvwS3njDTU/qOfwdd7j2g3fegdatvYb3wtoUjDEmSb9+8PbbbuGCu/3PwvP0024gdJ8+0KOH9/BeWFIwxhiA4cOhd2+4+WZ48UXv4d95B55/Hv7v/+CZZ7yH98aSgjHGfPON+7Zu2tT9lPc8emzCBDf2rWVLlxw8j33zypKCMSZvW7wYrr/edQEaO9atj+DRf//rxr7FxYVl7Jt3lhSMMXnX+vVu9FjJkq71t3Rpr+GXL3c9W6tUcTNlFC/uNXxYZPOcZYwxYbJjhxuctnu3W7ygSpX0jzkBSQuzFSjgxr5VqOA1fNhYUjDG5D0HD7rupitWhGVps1273AXI1q0wcyaceabX8GFlScEYk7ccOQKdO7tlNEeMgGbNvIY/dMg1USxZ4qZL8rwwW9hZUjDG5C29erkW3379XPdTj44cgdtuc52ZPvoImjf3Gj5LWEOzMSbvGDQIXn7Z9Q995BHv4Xv1gpEjoW9fuPVW7+GzhCUFY0zeMGaMW+Py2mtdcvA8WODVV12+uftueOwxr6GzlCUFY0zu98MP0LEjNG4Mn3wC+fN7DT96NPTsCTfcAK+9lr0Hp6XHkoIxJndbvtzNPFe1qhtaXLSo1/AzZriqoosvdjNleM43Wc6SgjEm9/rrr6MHC5Qv7zX84sWuNurss+GLL6BIEa/hI8J6Hxljcqfdu91kQ5s3u+6nngcL/PGHG/tWsqTLN2XKeA0fMWG9UhCRFiLym4isEpFjml5EpKqIzBCR/4nIYhG5OpzlMcbkEUkr2SxaBJ995iYe8mjbNncBsncvTJnifTB0RIUtKYhIfuBN4CqgDtBBROqk2O1J4FNVbQC0B94KV3mMMXmEKtx+u/v5PniwG1rs0f79roli9WoYP977YOiIC+eVQiNglaquUdVDwCigTYp9FCgV3D8F2BjG8hhj8oJnn4WhQ91KNt26eQ19+LDrxPTf/7pG5Usv9Ro+Wwhnm0IlYH3I43jgvBT7PAN8LSL3AsWBy8NYHmNMbjdkiEsKXbp4X8kmaaXOcePcmISbbvIaPtsI55VCaj11NcXjDsCHqloZuBr4WESOKZOI9BCReSIyb/PmzWEoqjEmx5syxVUbXXllWFayeekleOstNxD6/vu9hs5WwpkU4oHQ5pfKHFs91BX4FEBVfwSKAMf0GVPVd1U1TlXjKuSU+WeNMVln3jy48UaoX981LBcs6DX8Rx/B44+7qqOXXvIaOtsJZ1L4GaghItVFpBCuIXlCin3+AC4DEJHauKRglwLGmIxbs8Z1PT31VDctacmSXsN/9RV07QqXXx6WlTqznbC9PVVNBO4BpgLLcL2MlojIcyLSOtjtIaC7iCwCRgK3qWrKKiZjjEndli1usEBCgqs+Ou00r+GTLkDq1XNTJxUq5DV8thTWwWuqOhmYnGJbn5D7S4ELw1kGY0wuldQ3dN06mD4dzjnHa/jVq11v1goVXL4pVSr9Y3IDG9FsjMl5kvqGzpnj2hAuushr+E2bXHv1kSOu+sjzBUi2ZknBGJOzqLopsMeNc1OS3nCD1/B79rgmio0b4dtvoVYtr+GzPUsKxpic5eWX4Y034KGH3MABj5Jmx1iwwOWcxo29hs8RLCkYY3KOTz6BRx9139z/+Y/X0KrQo4drP3jnHddckRfl8s5VxphcY8YMtwDyJZfAsGHe+4b26QMffuj+7dHDa+gcxZKCMSb7++UXt3BBzZquXqdwYa/hBw+GF15wUyV5nh0jx7GkYIzJ3uLj3ViEEiVg8mTvCxeMH+/WVW7ZEt5+O2cvpemDtSkYY7KvnTtdQti1C2bPdktqevTDD9Chg1tuYfRot0BbXmenwBiTPR06BNdd59ZY/uorN6+RR8uWwTXXuAVyJk6E4sW9hs+xLCkYY7KfI0fc9NczZsDHH8Nll3kNv3GjWzmtUCGXb2yezX9YUjDGZD+PP+66n/btC506eQ2dVCO1bVtYlm7O8SwpGGOylzffhH794M474bFjlnY/KQcPwvXXw9KlbkLV2Fiv4XMFSwrGmOxj/Hi49143cuz11712BTpyxA1z+PZbN8yheXNvoXMV65JqjMkefvzRdQVq1AhGjoT8+b2Gf+QRGDXKLZJzyy1eQ+cqlhSMMZG3YoXrClS5Mnz5JRQr5jX8wIEwYIC7CHn0Ua+hcx1LCsaYyPr7b9cVKF++sHQFGjUKHnzQLZYzcKANTkuPtSkYYyJnzx5o1Qr++gtmzoSzzvIa/ttv4dZboUkT17PVc41UrmRJwRgTGYmJ0K6dm6f6iy9cW4JHixa5sW81a7r26yJFvIbPtSwpGGOynqrrcjp5spunulUrr+HXrXNjEUqVcjVSnqdLytWO26YgIjcF/1bPmuIYY/KE55+H99+HJ5/0Pk/11q2uiWLfPrc2QuXKXsPneuk1NPcO/h0T7oIYY/KIoUPh6addZf9zz3kNvX+/G+KwZg1MmABRUV7D5wnpVR9tE5EZQHURmZDySVXNo2sTGWMyZepU6N4drrgC3nvPa1egw4fh5pvdcIfRo13jsjlx6SWFq4FY4GPglfAXxxiTay1Y4PqFRkXB55+72eg8UYV77nENyq+9Bjfd5C10npNeUhiiqreIyHuqOitLSmSMyX3WrnWr2JQp4xqXS5XyGr5vX7d62qOPwn33eQ0JT64zAAAbRklEQVSd56TXptBQRM4AOopIGREpG3rLigIaY3K4bdtcy++BA64r0Omnew0/dKhrr+7UCV580WvoPCm9K4XBwFfAmcB8ILQCUIPtxhiTugMHXMvv77/DtGlQp47X8FOm/NNEMWSIGxRtTs5xT6GqDlLV2sAHqnqmqlYPuaWbEESkhYj8JiKrRCTVOXBFpK2ILBWRJSLySSbfhzEmuzl82P18/+EHN5zYc8vvzz+7JoroaBgzxmsTRZ6WocFrqnqniNQHLg42faeqi493jIjkB94ErgDigZ9FZIKqLg3Zpwau2+uFqrpdRE7NzJswxmQzqm7CoTFj4JVXoG1br+FXrXJNFKee6pooSpb0Gj5Py9DFlojcB4wATg1uI0Tk3nQOawSsUtU1qnoIGAW0SbFPd+BNVd0OoKqbTqTwxphsasAAGDQIHnjAJQePNm1yTRRHjrgerqed5jV8npfRaS66Aeep6l4AEekH/Ai8fpxjKgHrQx7HA+el2KdmEO8HID/wjKp+lTKQiPQAegBUrVo1g0U2xkTEqFHw8MOubucVvz3Z9+xxVwgbN7rJ7mrW9BrekPGpswU4HPL4MEc3Oqd1TEqa4nEBoAZwKdABeF9ESh9zkOq7qhqnqnEVbIVtY7KvmTOhc2e46CLXjuCx5TchwY0/+N//4NNPoXFjb6FNiIxeKQwFfhKRccHja4Eh6RwTD1QJeVwZ2JjKPnNUNQH4XUR+wyWJnzNYLmNMdrFkCVx7rZv++osvvE5Lqup6GX31lRsI7Xn+PBMiQ2lcVQcAXYBtwHagi6q+ms5hPwM1RKS6iBQC2gMpp8oYDzQFEJHyuOqkNRkvvjEmW9iwwVX0Fy3q+omW9TuM6amn4KOP4JlnoFs3r6FNChm6UhCRl4Ghqjooo4FVNVFE7gGm4toLPlDVJSLyHDBPVScEzzUXkaW4KqlHVHXrCb8LY0zk7NwJV18NO3bA7Nlwxhlew7/9Nvz73+5KoU8fr6FNKkQ1ZTV/KjuJdMNdKRTAVSWNVNWdYS5bquLi4nTevHmReGljTEqHDrmEMGsWTJoEzZt7DT9uHNxwg6suGjsWCtgKMJkmIvNVNS69/TJaffS+ql4I3ApUAxaLyCci0vTkimmMybFUoWtX+OYbtzaC54Tw/ffQoQOcd57r0GQJIWtkuGtAMBjtnOC2BVgEPCgio8JUNmNMdvbkkzB8uFswp3Nnr6GXLoVrrnE1UV9+CcWKeQ1vjiOjbQoDgNbAN0BfVZ0bPNUv6DFkjMlLBg92U5P26AFPPOE1dFKbdZEirrdR+fJew5t0ZPSC7FfgSVXdl8pzflfbNsZkbxMmwN13u4r+N9/0ulDOjh1ubeXt2+G776C6LQSc5TI699EHwdTZUUCRkO3fRarB2RgTAXPmQPv20LCh94r+gwfhuutg2TI3n1GDBt5CmxOQ0eqjbsD9uAFoC4HGuGkumoWvaMaYbGXlSlfRX7EiTJwIxYt7C33kiFuyeeZMNxD6iiu8hTYnKKMNzfcD5wLrVLUp0ADYHLZSGWOyl02bXL0OuIr+U/1OaPzww27qiv/8x822bSIno9d+B1T1gIggIoVVdbmI1ApryYwx2cPeva79YONGmDEDatTwGn7AABg4EO6/3yUHE1kZTQrxwUR144FpIrKdY+cxMsbkNomJrg1h/nw3kuy8lBMdn5yRI+Ghh9xEdwMGeG2zNpmU0Ybm64K7z4jIDOAU3DKdxpjcShXuuce1H7z1lltW06NvvnHDGy65BIYNs6U0s4vjJgURKQLcAZwN/AIMUdVZWVEwY0yE9e0L77wDjz0Gd97pNfTCha6nUa1aMH681wlVzUlKLzd/BMThEsJVgN8VM4wx2dNHH7kRy506ueTg0dq1brqkU05xE6qWPmYFFRNJ6VUf1VHVegAiMgSYm87+xpicbto0Nz/1ZZfBkCFeK/q3bnWjlffvd3MbVa7sLbTxJL2kkJB0J5gKO8zFMcZE1MKFcP31UKcOjBkDhQp5C71/vxvmsHatyzt163oLbTxKLynUF5FdwX0BigaPBVBVLRXW0hljss66dW4sQunSbkjxKad4C52Y6GY8nTMHPvsMLr7YW2jj2XGTgqrmz6qCGGMiaNs2lxD274cffoBKlbyFVoV773UrdL7+ulsfwWRfNkO5MXndgQNubeXVq2HqVO/1Ov/+t5tU9bHHXA9Xk71ZUjAmL0uadGj2bDeS7NJLvYYfOtStr3zrrd47MZkwseEixuRljzziKvlfftmNXPZo8mS3rvKVV7qF2ayfSs5gScGYvOrVV93cEvfdBw8+6DX03Llu6oqYGPj8cyhY0Gt4E0aWFIzJiz77zCWC66/3PunQypXQsiWcdhpMmgQlSngLbbKAJQVj8prZs+GWW+CCC9way/n9dTL8+283OA3cDNv/+pe30CaLWEOzMXnJsmVuYrtq1Vwf0aJFvYXevdtNX/HXX2GZYdtkEUsKxuQVGze6n/FFirif8eXKeQudkODaEBYtcrmmka3cnmNZUjAmL9i1y1X0b90Ks2a5KwVPVN1USVOnul5GLVt6C20iIKxtCiLSQkR+E5FVIvLYcfa7UURUROLCWR5j8qSEBLjxRvjlF9cVqGFDr+GfeMKth/Dss9C1q9fQJgLClhREJD/wJm7K7TpABxGpk8p+JYH7gJ/CVRZj8ixVN1hg2jR4771/WoE9eeMNePFF6NHDDVIzOV84rxQaAatUdY2qHgJGAW1S2e954D/AgTCWxZi8qU8ftzbCs89Cly5eQ48d64Y4tG4Nb75pg9Nyi3AmhUrA+pDH8cG2ZCLSAKiiqhOPF0hEeojIPBGZt3nzZv8lNSY3eucdeOEFV+Hv+Wf87Nlw883QuLGbHaOAtU7mGuFMCqn9btDkJ0XyAQOBh9ILpKrvqmqcqsZVqFDBYxGNyaUmToS77nJ9RN9+2+vP+KVL/+nV+uWXUKyYt9AmGwhnUogHqoQ8rgxsDHlcEogCZorIWqAxMMEam405SXPnQrt20KABjB7t9Wd8fHzYerWabCKcF30/AzVEpDqwAWgP3Jz0pKruBMonPRaRmcDDqjovjGUyJndbtQpatXJDiT3PMbFzp7vw2LEDvvvOa69Wk42E7UpBVROBe4CpwDLgU1VdIiLPiUjrcL2uMXnW5s1uoZwjR2DKFK9zTBw86JZcWL7cNTDHxHgLbbKZsDYPqepkYHKKbX3S2PfScJbFmFxt3z63AHJ8PHzzDdSq5S100pILM2fCiBFw+eXeQptsyPoMGJPTHT7sFkCeOxfGjHET3Xmi6iZT/fRT6N/f9TgyuZslBWNysqQFkCdMcAsgX3ed1/CvvAKvvQb33w8PpdtP0OQGNnW2MTlZv36uy+mjj3pfAHnECLcw2003eV9ywWRjlhSMyamGD4fevV3V0Ysveg09fbobAH3JJW5eo3z2TZFn2J/amJzom2/g//4PLr0Uhg71+q29cKFbkO2cc2D8eDcmweQdlhSMyWkWL3bf2rVqwbhxULiwt9C//+56tZYu7Xq1li7tLbTJIayh2ZicZP16961dsqT3b+0tW9xo5QMH3IVIpUrpH2NyH0sKxuQUO3a4hLBnD3z/PVSu7C30vn1uPqN161x7Qp1jJrk3eYUlBWNygqQhxStWuCXO6tXzFjox0bVVz5nj1uC56CJvoU0OZEnBmOzuyBG47Ta3jOaIEdC0qbfQqnD33W6Yw5tvuqYKk7dZQ7Mx2V2vXjBqlBuT4HlI8fPPw7vvwuOPu5m2jbGkYEx2NmgQvPyyG5j2yCNeQw8ZAk8/DZ07u7V4jAFLCsZkX2PGwAMPuLaEV1/1OqR40iS4/Xa48kq3dLONVjZJLCkYkx398AN07OjWu/zkE8if31von35yU1fExLiG5YIFvYU2uYAlBWOym+XL3TTYVau6FuCiRb2FXrnSrcFTsaL3NXhMLmFJwZjs5K+/3AiyggXdepfly6d/zAmEvvJKd/+rr7yuwWNyEeuSakx2sXs3tGzphhbPnAlnnuk99N9/w4wZUKOGt9Aml7GkYEx2kJDgKvoXLXJVRnFx3kIfOgQ33PBP6EaNvIU2uZAlBWMiTdV1BZo61XUFuvpqr6G7doVp0+CDD7yGNrmUtSkYE2nPPuumv+7TB7p18xq6d2+37MLzz7v1EYxJjyUFYyJpyBCXFLp0gWee8Rr69dfdIOg77oAnnvAa2uRilhSMiZQpU/4ZQfbOO15HkI0Z49ZVbtMG3njDBqeZjLOkYEwkzJvnGpajo+Gzz7yOIJs92417O/98GDnS67g3kwdYUjAmq/3+uxtBVr68G0FWsqS30EuWuHURqlf3Pu7N5BHW+8iYrLR1qxucduiQG4tQsaK30PHxLnTRom5wWrly3kKbPCSsVwoi0kJEfhORVSLyWCrPPygiS0VksYh8IyJnhLM8xkTU/v1u+op169zP+HPO8RY6aVG2XbtcU8UZ9j/JZFLYkoKI5AfeBK4C6gAdRCTlIn//A+JUNRr4HPhPuMpjTEQdPuwq+ufMcQvleFze7MAB16D8228wbhzUr+8ttMmDwnml0AhYpaprVPUQMApoE7qDqs5Q1X3BwzmAv0VnjckuVN0U2OPGuSmwb7jBW+jDh+GWW+C772DYMGjWzFtok0eFMylUAtaHPI4PtqWlKzAltSdEpIeIzBOReZs3b/ZYRGOywMsvu36hDz0E993nLawq9Ozppr9+5RVo395baJOHhTMppNYzWlPdUaQTEAf0T+15VX1XVeNUNa5ChQoei2hMmI0cCY8+Cu3awX/81o727+8GqPXsCQ8+6DW0ycPC2fsoHqgS8rgysDHlTiJyOfAEcImqHgxjeYzJWjNmuLUuL7kEPvoI8vn7DTZ8uFu6uV07dyFijC/hvFL4GaghItVFpBDQHpgQuoOINADeAVqr6qYwlsWYrPXLL3DddVCzpmtLKFzYW+hp09ysGE2bes81xoQvKahqInAPMBVYBnyqqktE5DkRaR3s1h8oAXwmIgtFZEIa4YzJOeLj3XSkxYvD5MlQpoy30AsWwPXXQ5063nONMUCYB6+p6mRgcoptfULuXx7O1zcmy+3c6RLCzp1uvomqVb2F/v13F7pMGTcW4ZRTvIU2JpmNaDbGl0OHXJXRsmXuW9vjgIEtW9y8eYcOuaaK00/3FtqYo1hSMMaHI0dcRf+MGW7AwOX+LoL37nVTJa1fD9OnQ+3a3kIbcwxLCsb48Pjj8Mkn0LevG03mSWKiG3/w889uOuwLL/QW2phUWVIw5mS99dY/q9k8dswUX5mmCnfeCRMnupe49lpvoY1Jk3VmM+ZkjB8P99zjJrp7/XWvq9k89xy8/75bNe3OO72FNea4LCkYk1k//ggdOkCjRjBqFBTwd+H9/vtudc7bbnPrKxuTVSwpGJMZK1a4q4PKleHLL6FYMW+hJ050q3RedRW8+64tpWmyliUFY07U33+71Wzy5XOr2Xicj2vOHGjbFmJj4dNPva7SaUyGWEOzMSdizx5o2dIlhhkz4KyzvIVescJ1PT39dLdKZ4kS3kIbk2GWFIzJqMRE9zP+f/+DL75wbQme/PWXG5yWdPFx6qneQhtzQiwpGJMRSf1Dp0yBwYPdT3pPdu9201ds3uyWbT77bG+hjTlhlhSMyYjnn/+nf+jtt3sLe+iQW4ht8WLXXh0X5y20MZliScGY9AwdCk8/Dbfe6rV/6JEj0LWrmwp76FDX28iYSLPeR8Ycz9Sp0L07XHEFvPee1/6hvXu7xXJeeMGNRzAmO7CkYExaFixwdTtRUW4h5EKFvIUeNMitznnXXW7aJGOyC0sKxqRm7VrX+luunFsop1Qpb6E//xweeMDNZTRokA1OM9mLtSkYk9K2bW5w2sGD8O23XhcvmDULOnaECy5wk6rmz+8ttDFeWFIwJtT+/dC6tVvmbNo0t+6lJ7/8Am3auPFuEyZA0aLeQhvjjSUFY5IcPgydOsEPP8Do0dCkibfQ69e73kXFi7vBaWXLegttjFeWFIwBNzjtwQdh7FgYMMCNXPZk+3ZXG7V7t/dlm43xzpKCMeASwaBBrgW4Z09vYQ8ccFVGK1e63q3R0d5CGxMWlhSMGTUKHn4YbrwRXnnFW9jDh12j8uzZMHIkNG3qLbQxYWNdUk3eNmsWdO4MF10EH3/sZqTzQNVddCTVRrVv7yWsMWFnScHkXUuWuMECZ57pZj0tUsRb6H794I03XDOFx9ooY8LOkoLJm37/3XUHKlLEzXzqsTvQsGFuCosOHaB/f29hjckSYU0KItJCRH4TkVUi8lgqzxcWkdHB8z+JSLVwlscYwK1m06SJWzBn8mSoVs1LWFXXPNG1KzRr5ia581QbZUyWCdtHVkTyA28CVwF1gA4iknIkUFdgu6qeDQwE+oWrPMYA8OuvLiEcPOhWTmvQwEvY+fNdQ3KHDhAT49oSChf2EtqYLBXO3keNgFWqugZAREYBbYClIfu0AZ4J7n8OvCEioqrquzDz/j2V0s8/6DusyWEqJvzB3nwl6XrGDNbcXNtLzCNHYPlyt1TzW2+5SVULWL8+k0OF86NbCVgf8jgeOC+tfVQ1UUR2AuWALaE7iUgPoAdA1UyO/ClcoRSby/ubssDkTPEF4hhT+ymKlDgbn5+Gtm1do/Ipp3gMakwEhDMppDb3Y8orgIzsg6q+C7wLEBcXl6mriHo9zocen2XmUJPLXBrpAhiTjYWzGSweqBLyuDKwMa19RKQAcAqwLYxlMsYYcxzhTAo/AzVEpLqIFALaAxNS7DMB6BzcvxH4NhztCcYYYzImbNVHQRvBPcBUID/wgaouEZHngHmqOgEYAnwsIqtwVwg27tMYYyIorH0kVHUyMDnFtj4h9w8AN4WzDMYYYzLOhtYYY4xJZknBGGNMMksKxhhjkllSMMYYk0xyWg9QEdkMrMvk4eVJMVo6m8iu5YLsWzYr14mxcp2Y3FiuM1S1Qno75bikcDJEZJ6qxkW6HCll13JB9i2blevEWLlOTF4ul1UfGWOMSWZJwRhjTLK8lhTejXQB0pBdywXZt2xWrhNj5ToxebZceapNwRhjzPHltSsFY4wxx2FJwRhjTLIcmxRE5AMR2SQiv6by3MMioiJSPo1jO4vIyuDWOWR7QxH5RURWicggEUltEaCwlEtEYkTkRxFZIiKLRaRdyHMfisjvIrIwuMVkVbmC5w+HvPaEkO3VReSn4DyODqZIz5JyiUjTkDItFJEDInJt8FxYzpeIPCMiG0LiXp3GsS1E5Lfgc/RYyPawnK+MlEtEqojIDBFZFnzG7j/R9xWusgX7rQ3+7y0UkXkh28uKyLTgnE0TkTJZUSYRqZXi87VLRB44kfeUmbIF2+8NPj9LROQ/aRwbts8Yqpojb0ATIBb4NcX2KrjputcB5VM5riywJvi3THC/TPDcXOB83IpwU4CrsrBcNYEawf3TgT+B0sHjD4EbI3G+gn32pLH9U6B9cH8wcGdWlivF33QbUCyc5wu3nvjD6RyXH1gNnAkUAhYBdcJ5vjJYropAbHC/JLAipFzpHh/OsgX7rU3j/8V/gMeC+48B/bKqTCn+pn/hBn+F+3w1BaYDhYPHp2b1ZyzHXimo6nekvkrbQOBRUlnWM3AlME1Vt6nqdmAa0EJEKgKlVPVHdWd0GHBtVpVLVVeo6srg/kZgE5Du6MNwlystIiJAM+DzYNNHZOH5SuFGYIqq7jvR189EudLTCFilqmtU9RAwCmiTBecrveP+VNUFwf3dwDLcGunenMQ5O542uHMFmThnnsp0GbBaVTM7k0Kq0ijbncBLqnow2GdTKoeG9TOWY5NCakSkNbBBVRcdZ7dKwPqQx/HBtkrB/ZTbs6pcofs3wv0CWB2y+d/iqpUGikjhLC5XERGZJyJzkqpogHLADlVNDB5H7HzhFmcamWKb9/MVuCeI+0EaVRlpfb7Cdr4yWK5kIlINaAD8lJnjw1Q2Bb4Wkfki0iNk+79U9U9wiQ04NQvLlCS1z1e4zldN4OKgCmiWiJybyj5h/YzlmqQgIsWAJ4A+6e2ayjY9zvasKlfS/hWBj4Euqnok2NwbOAc4F1dV0iuLy1VV3dD6m4FXReQsstf5qoerakri/XwF3gbOAmJw1XuvpFakVLaF7fN1AuUCQERKAGOAB1R114keH8ayXaiqscBVwN0i0sRjGTJbJoI6+dbAZ5k5PhMK4Kq1GwOPAJ8GVwBHFSuV47x9xnJNUsD9kaoDi0RkLVAZWCAip6XYLx5XX52kMrAx2F45le1ZVS5EpBQwCXhSVeckbQ8u/TW4pByKu3zMsnIF1Vmo6hpgJu5X5hagtIgkrd6X5ecr0BYYp6oJIeUNx/lCVf9W1cNBsn4vjbhpfb7Cdb4yWi5EpCAuIYxQ1bEnenw4yxbyGdsEjAvZ7+8g8Sf9AEitOiUsZQpcBSxQ1b8zefyJigfGBp/fucAR3CR4KfcJ22cs1yQFVf1FVU9V1WqqWg134mJV9a8Uu04FmotImeCyrzkwNbg03S0ijYPMfCvwRVaVK/hFMg4YpqqfpXgu6T+F4OoIj+mpE8ZylUmqfhHXC+hCYGnQ7jIDV58P0JksPF8hOpDi0j4c5ys0buC6NOL+DNQIeoEUwlU9TAjX+cpouYJzMQRYpqoDTvT4MJetuIiUTLqP+z+ZtN8E3LkCT+fsBN9vmp+vDB5/osbj2gUQkZq4auSUs6KG9zN2oi3T2eWG+0P9CSTgvji6pnh+LUFvBiAOeD/kuf8DVgW3LiHb43B/4NXAGwQjvrOiXECn4JiFIbeY4LlvgV+Csg0HSmRhuS4IXntR8G/XkGPOxPXYWoW7vC6cxX/HasAGIF+KY8JyvnDVer8Ai3FfVhWDfU8HJoccezWud89q4Ilwn6+MlAu4CFeVsDjk83V18Fyqx2dh2c4MPl+LgCUpzlk54BtgZfBv2Sz8OxYDtgKnpIgZzvNVKPjM/gosAJpl9WfMprkwxhiTLNdUHxljjDl5lhSMMcYks6RgjDEmmSUFY4wxySwpGGOMSWZJwRgPROQJ+WeG24Uicl6ky2RMZhRIfxdjzPGIyPlAK9wgu4PBIL8Tn7LYmGzAkoIxJ68isEX/mdky5QhUY3IMG7xmzEkKJpn7HjcCdjowWlVnRbZUxmSOtSkYc5JUdQ/QEOgBbAZGi8htES2UMZlkVwrGeCYiNwKdVfWaSJfFmBNlVwrGnCRx6/nWCNkUg1tG1JgcxxqajTl5JYDXRaQ0kIibobLH8Q8xJnuy6iNjjDHJrPrIGGNMMksKxhhjkllSMMYYk8ySgjHGmGSWFIwxxiSzpGCMMSaZJQVjjDHJ/h+9PSeNhMmXcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "S = np.linspace(14,16,100)\n",
    "K = 15            # strike\n",
    "cost_b = 0.01     # transaction cost \n",
    "\n",
    "plt.plot(S, np.maximum(S-K,0), color='blue',label=\"Zero cost Payoff\")\n",
    "plt.plot(S, np.maximum(S*(1+cost_b)-K,0), color='red',label=\"Transaction costs Payoff\")\n",
    "plt.xlabel(\"S\")\n",
    "plt.ylabel(\"Payoff\")\n",
    "plt.title(\"Payoff comparison\")\n",
    "plt.legend(loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec2'></a>\n",
    "## Singular control problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec2.1'></a>\n",
    "### Maximization problem\n",
    "\n",
    "The **value function** of the maximization problem for $j=0,w,b$ (corresponding to the three portfolios: no option, writer, buyer) is defined as:\n",
    "\n",
    "\\begin{equation}\n",
    "V^j(t,b,y,s) = \\sup_{L,M} \\;  \\mathbb{E}\\biggl[ \\; \\mathcal{U}( \\mathcal{W}^{j}_T ) \\; \\bigg| \\; B_{t} = b, Y_{t} = y, S_{t} = s \\biggr],             \n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathcal{U}: \\mathbb{R} \\to \\mathbb{R}$ is a concave increasing **utility function**. The **exponential utility** is what we are looking for:\n",
    "\n",
    "\\begin{equation}\n",
    " \\mathcal{U}(x) := 1- e^{-\\gamma x}   \\quad \\quad \\gamma >0 .\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec2.2'></a>\n",
    "### Indifference pricing\n",
    "\n",
    "The writer (buyer) option price is defined as the amount of cash to add (subtract) to the bank account, \n",
    "such that the maximal expected utility of wealth of the writer (buyer) is the same he could get with \n",
    "the zero-option portfolio.\n",
    "\n",
    "* The **writer price** is the value $p^w>0$ such that \n",
    " \\begin{equation}\n",
    "  V^0(t,b,y,s) = V^w(t,b+p^w,y,s),\n",
    " \\end{equation}\n",
    " \n",
    "* The **buyer price** is the value $p^b>0$ such that\n",
    " \\begin{equation}\n",
    "  V^0(t,b,y,s) = V^b(t,b-p^b,y,s).\n",
    " \\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec3'></a>\n",
    "## Variable reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the properties of the exponential utility, it is possible to remove $\\mathbf{B}$ from the state variables.\n",
    "\n",
    "$$ V^j(t,b,y,s) = \\sup_{L,M} \\; \\mathbb{E}_{t,b,y,s}\\biggl[  1- e^{-\\gamma \\mathcal{W}^j(T) } \\biggr]   \n",
    "\t     = 1- e^{-\\gamma \\frac{b}{\\delta(t,T)}} Q^j(t,y,s), $$\n",
    "         \n",
    "where $\\delta(t,T) = e^{-r(T-t)}$. (for the full calculations, check the paper. Equations 4.21 -4.25).\n",
    "\n",
    "<a id='sec3.1'></a>\n",
    "### Minimization problem\n",
    "\n",
    "\\begin{equation}\n",
    "Q^j(t,y,s) = \\inf_{L,M} \\; \\mathbb{E}_{t,y,s}\\biggl[ \\;\n",
    "\t     e^{-\\gamma \\bigl[ -\\int_{t}^T (1+\\theta_b) \\frac{S_u}{\\delta(u,T)} dL_u +\n",
    "\t     \\int_{t}^T (1-\\theta_s) \\frac{S_u}{\\delta(u,T)} dM_u \\bigr] } \\; H^j(Y_T,S_T) \\bigg]  \n",
    "\\end{equation}\n",
    "\n",
    "The exponential term inside the expectation can be considered as a discount factor, and the second term is the terminal payoff:\n",
    " - No option:\n",
    " \n",
    " $$ H^0(y,s) = e^{-\\gamma \\, c(y,s)}. $$\n",
    " \n",
    " - Writer:\n",
    " \n",
    " $$ H^w(y,s) = e^{-\\gamma \\bigl[ c(y,s)\\mathbb{1}_{\\{s(1+\\theta_b) \\leq K\\}} + \n",
    " \\bigl( c( y-1,s) + K \\bigr) \\mathbb{1}_{\\{s(1+\\theta_b)>K\\}} \\bigr] }.$$\n",
    " \n",
    " - Buyer:\n",
    "\n",
    "$$  H^b(y,s) = e^{-\\gamma \\bigl[ c(y,s)\\mathbb{1}_{\\{s(1+\\theta_b) \\leq K\\}} + \n",
    " \\bigl( c( y+1,s) - K \\bigr) \\mathbb{1}_{\\{s(1+\\theta_b)>K\\}} \\bigr] }.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec3.2'></a>\n",
    "### Portfolio dynamics (2 state variables)\n",
    "\n",
    "In order to simplify the numerical computations,let us pass to the log-variable $X_t = \\log(S_t)$.\n",
    "\n",
    "The resulting portfolio dynamics is:\n",
    "\n",
    "\\begin{equation}\n",
    " \\begin{cases}\n",
    " dY_t &=  dL_t - dM_t \\\\\n",
    " dX_t &= \\biggl( \\mu - \\frac{1}{2} \\sigma^2 \\biggr) dt + \\sigma dW_t.\n",
    "\\end{cases}\n",
    "\\end{equation} \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec3.3'></a>\n",
    "### HJB variational inequality\n",
    "\n",
    "The Hamilton Jacobi Bellman equation associated to the minimization problem is:\n",
    "\n",
    "$$\n",
    " \\min \\; \\biggl\\{ \\; \\frac{\\partial Q^j}{\\partial t} + (\\mu-\\frac{1}{2}\\sigma^2) \\frac{\\partial Q^j}{\\partial x}\n",
    "+ \\frac{1}{2}\\sigma^2 \\frac{\\partial^2 Q^j}{\\partial x^2}  ,\n",
    " \\; \\frac{\\partial Q^j}{\\partial y} +(1+\\theta_b) e^x \\frac{\\gamma}{\\delta(t,T)}Q^j \\; , \n",
    "\\; -\\biggl( \\frac{\\partial Q^j}{\\partial y}+(1-\\theta_s)e^x \\frac{\\gamma}{\\delta(t,T)} Q^j \n",
    "\\biggr) \\biggr\\} = 0. \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec3.4'></a>\n",
    "### Indifference price (explicit form)\n",
    "\n",
    "Using again the explicit form of the utility function, we obtain formulas for the option prices:\n",
    "\n",
    "$$\n",
    " p^w(t_0,y,x) = \\frac{\\delta(t_0,T)}{\\gamma} \\log \\biggl( \\frac{Q^w(t_0,y,e^x)}{Q^0(t_0,y,e^x)} \\biggr), $$\n",
    " \n",
    "$$ p^b(t_0,y,x) = \\frac{\\delta(t_0,T)}{\\gamma} \\log \\biggl( \\frac{Q^0(t_0,y,e^x)}{Q^b(t_0,y,e^x)} \\biggr).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec4'></a>\n",
    "# Numerical Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec4.1'></a>\n",
    "###  Discrete SDE\n",
    "\n",
    "As usual, we introduced the time step $\\Delta t = \\frac{T}{N}$, where we assumed $t_0 = 0$ and $N \\in \\mathbb{N}$. \n",
    "The time $t_n = n \\Delta t$, for $n \\in \\{0,1,2, ..., N\\}$. \n",
    "\n",
    "The space discretization has 2 dimensions:\n",
    "- The space step $h_x$ is defined as $h_x = \\sigma \\sqrt{\\Delta t}$.\n",
    "- The space step is $h_y$. In this computations we choose $h_x = h_y$.\n",
    "\n",
    "The discretized version of the Stochastic Differential equation is: \n",
    "\n",
    "$$\n",
    " \\begin{cases}\n",
    " \\Delta Y_n &= \\; \\Delta L_n - \\Delta M_n \\\\\n",
    " \\Delta X_n &= \\; (\\mu - \\frac{1}{2} \\sigma^2 )  \\Delta t + \\sigma \\Delta W_n\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Both $\\Delta L_n$ and $\\Delta M_n$ at each time $t_n$ can assume values in $\\{0,h_y\\}$. They cannot be different from zero at the same time (It is quite strange to buy and sell at the same time, right?)\n",
    "\n",
    "The variable $\\Delta W_n$ has $\\frac{1}{2}$ probability of being equal to $h_x$ and $\\frac{1}{2}$ probability of being equal to $-h_x$.\n",
    "\n",
    "The variation $\\Delta X_n$ is $\\pm h_x$ plus the drift component. We obtain a recombining **binomial tree**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Binomial tree with drift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.7080502]\n",
      "[2.61744646 2.82157061]\n",
      "[2.52684272 2.73096687 2.93509101]\n",
      "[2.43623898 2.64036313 2.84448727 3.04861142]\n",
      "[2.34563524 2.54975939 2.75388353 2.95800768 3.16213182]\n",
      "[2.2550315  2.45915565 2.6632798  2.86740394 3.07152809 3.27565223]\n"
     ]
    }
   ],
   "source": [
    "N = 6; dt=1/N; S0 = 15; x0 = np.log(S0) \n",
    "mu = 0.1; sig = 0.25; h_x = sig * np.sqrt(dt)\n",
    "\n",
    "for n in range(N):\n",
    "    x = np.array( [x0 + (mu-0.5*sig**2)*dt*n + (2*i-n)*h_x for i in range(n+1) ] )\n",
    "    print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec4.2'></a>\n",
    "###  Algorithm\n",
    "\n",
    "Using the Dynamic Programming Principle (DPP) on the minimization problem we obtain a recursive algorithm on the nodes of the grid.\n",
    "\n",
    "$$ \\begin{aligned}\n",
    " Q^{j}(t_n,Y_n,X_n) = \\min  \n",
    " & \\; \\biggl\\{ \\mathbb{E}_n \\biggl[ Q \\bigl( t_{n+1}, Y_n, X_n + \\Delta X_n \\bigr) \\biggr], \\\\ \\nonumber\n",
    " & \\; \\exp \\biggl(\\frac{\\gamma}{\\delta(t_n,T)} (1+\\theta_b) e^{X_n} \\Delta L_n \\biggr) \n",
    "  \\mathbb{E}_n \\biggl[ Q^{j} \\bigl( t_{n+1}, Y_n+\\Delta L_n, X_n + \\Delta X_n \\bigr) \\biggr], \\\\ \\nonumber\n",
    " & \\; \\exp \\biggl(\\frac{-\\gamma}{\\delta(t_n,T)} (1-\\theta_s) e^{X_n} \\Delta M_n \\biggr)\n",
    "  \\mathbb{E}_n \\biggl[ Q^{j} \\bigl( t_{n+1}, Y_n-\\Delta M_n, X_n + \\Delta X_n \\bigr) \\biggr]\n",
    " \\biggr\\}.\n",
    "\\end{aligned} $$\n",
    "\n",
    "#### How to solve this problem?\n",
    "- Create a **binomial tree** with N time steps. \n",
    "- Create a \"shares vector\" **y** with dimension M. \n",
    "- Initialize a two dimensional grid of size $(N+1)\\times M$, to be filled with the values of the terminal conditions $H^j(y,e^x)$ for $j=0,w,b$ (see [Minimization problem](#sec3.1))\n",
    "- Create a backward loop over time and find $Q^j(0,0,X_0)$\n",
    "\n",
    "#### Computational complexity?  Well... Quite high. \n",
    "- A binomial tree with N time steps has $\\sum_{n=0}^N (n+1) = \\frac{N(N+1)}{2} + (N+1) = (N+1)(\\frac{N}{2}+1)$ nodes. \n",
    "The loop over all the nodes is $\\mathcal{O}(N^2)$.\n",
    "- If we assume $M \\sim N$, the loop over all the values of **y** has another $\\mathcal{O}(N)$.\n",
    "- The minimum search for every point in **y**, produces another $\\mathcal{O}(N)$ operations.\n",
    "\n",
    "Therefore the total computational complexity is of $\\mathcal{O}(N^4)$.\n",
    "\n",
    "\n",
    "#### For space reasons, I will not expose the code here in the notebook. The interested reader can peek the (clear and commented) code inside the class TC_pricer.  [code](./functions/TC_pricer.py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec5'></a>\n",
    "# Numerical computations\n",
    "\n",
    "Let us import the classes we need:\n",
    " - **Option_param**: creates the option object\n",
    " - **Diffusion_process**: creates the process object\n",
    " - **TC_pricer**: creates the pricer for this model\n",
    " - **BS_pricer**: creates the Black-Scholes pricer, useful for making comparisons.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functions.Parameters import Option_param  \n",
    "from functions.Processes import Diffusion_process\n",
    "from functions.TC_pricer import TC_pricer\n",
    "from functions.BS_pricer import BS_pricer\n",
    "\n",
    "# Creates the object with the parameters of the option\n",
    "opt_param = Option_param(S0=15, K=15, T=1, exercise=\"European\", payoff=\"call\" )\n",
    "\n",
    "# Creates the object with the parameters of the process\n",
    "diff_param = Diffusion_process(r=0.1, sig=0.25, mu=0.1)\n",
    "\n",
    "# Creates the object of the Transaction Costs pricer\n",
    "TC = TC_pricer(opt_param, diff_param, cost_b=0, cost_s=0, gamma=0.0001)\n",
    "# Creates the object of the Black-Scholes pricer\n",
    "BS = BS_pricer(opt_param, diff_param)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We expect that if the transaction costs are **zero**, and the risk aversion coefficient $\\gamma \\to 0$ (i.e. the investor becomes risk neutral), the price should **converge** to the **Black-Scholes price**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Zero TC price:  2.246375063664713\n",
      "Black Scholes price: 2.246368616746695\n",
      "Difference: 6.446918018099268e-06\n"
     ]
    }
   ],
   "source": [
    "tc = TC.price(N=2000)\n",
    "bs = BS.closed_formula()\n",
    "\n",
    "print(\"Zero TC price: \", tc )\n",
    "print(\"Black Scholes price:\", bs )\n",
    "print(\"Difference:\", np.abs(tc-bs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Wait a second!!! WE CAN DO BETTER!\n",
    "\n",
    "##### Let us analyze the the writer and buyer prices, for different initial stock values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = list(range(5,21))\n",
    "price_0 = []; price_w = []; price_b = []\n",
    "\n",
    "for s in S:\n",
    "    TC.S0 = s\n",
    "    TC.cost_b=0; TC.cost_s=0\n",
    "    price_0.append(TC.price(N=400))   # zero costs\n",
    "    TC.cost_b=0.05; TC.cost_s=0.05\n",
    "    price_w.append(TC.price(N=400, TYPE=\"writer\")) \n",
    "    price_b.append(TC.price(N=400, TYPE=\"buyer\"))\n",
    "TC.cost_b=0; TC.cost_s=0  # set to 0 for future computations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd8TtcfwPHPyUASe1NitGpFRIgVe2vtPWqPGqWoVlVbrVFbtVU1iigpsUvUCkKIESP2Kk2CBImSyJB5fn/cJL8gIeO5TyLO+/V6XvLc5977PXnwfc5z7rnfI6SUKIqiKG8Pk8xugKIoimJcKvEriqK8ZVTiVxRFecuoxK8oivKWUYlfURTlLaMSv6IoyltGJX5Fd0IIJyHEjAyeo4kQ4q6h2pRCjIFCiKPpPPY7IcQ6Q7cpHe24LIRoktntULI2lfiVDBNC+AghIoQQoUKIx0KIXUKI0pnYnhxCiAVCiLvxbfpXCPFjZrXHmKSUVaWU7pndDiVrU4lfMZT2UsrcQAngAfBLJrZlMlALqA3kAZoC5zKxPboTQphldhuUN4dK/IpBSSmfAZuBKsm9LoQoIIRwFUIExn87cBVClEryekEhxGohhH/869tTOM9YIcSVpMcm4QBsk1L6S42PlPKPJMeWFkJsjW/DIyHE4hfOPT8+9r9CiLZJtpcUQuwQQvwnhPhHCDEspfdBCFFXCOEphHgihDifdPglfkjpthDiaXyMvimc4zshxGYhhEv8vmeFENWTvO4jhJgkhLgAhAkhzOK3tYh/3VQI8ZUQ4lb88WcSvokJISoJIfbH/y7XhRA9UvpdlOxHJX7FoIQQlkBP4EQKu5gAq4EygDUQASRNvGsBS6AqUBR4aYhGCPENMBBoLKVMbtz/BDBBCDFKCFFNCCGSHGsKuAK+QFngHWBDkmPrANeBwsBcYGWS49cDd4GSQDfgByFE82Ta9w6wC5gBFAQmAluEEEWEEFbAz0BbKWUeoD7gndwbFa8jsCn+PH8C24UQ5kle7w18COSXUsa8cOyE+Nc/APICg4Hw+Dbsjz9f0fh9lgghqr6iHUp2IqVUD/XI0APwAUKBJ0AM4A9US/K6EzAjhWPtgMfxP5cA4oACyezXBLgHLASOAvle0R5TYDRwDIiMb8+A+NfqAYGAWTLHDQT+SfLcEpBAcaA0EAvkSfL6LMAp/ufvgHXxP08C1r5w7r3AAMAq/n3qCli85n39DjiR5LkJEAA0TPK+D07m76JF/M/XgY7JnLcn4PHCtmXA1Mz+t6QexnmoHr9iKJ2klPmBnMAnwGEhRPEXdxJCWAohlgkhfIUQIcARIH98T7w08J+U8nEKMfIDw4FZUsrglBoipYyVUv4qpXSMP2YmsEoIUTk+hq98uXec4H6S84TH/5gbrZf/n5TyaZJ9fdG+MbyoDNA9fpjniRDiCdAAKCGlDENLvCOAgPgL4ZVS+l2AO0naE8f/v3G89HoySgO3UmhfnRfa1xftA055C6jErxhUfNLditY7bpDMLp8BFYE6Usq8QKP47QItiRUUQuRP4fSPgXbAaiGEYyrbEyGl/DX+2CrxMazTcTHUP75teZJss0b7FvKiO2g9/vxJHlZSytnxbdorpWyJ9g3nGrDiFXETZ0cJIUyAUvFtSfwVX3HsHeDdFLYffqF9uaWUI19xLiUbUYlfMSih6QgUAK4ms0setHH9J0KIgsDUhBeklAHAbrTx5gJCCHMhRKOkB0ttqmJfYJsQok4KbRgntHn/FvEXPAfExz0HnEIbLpkthLASQuRKzYeIlPIO4AnMij/GFhgCOCez+zqgvRCidfwF1lzx7SklhCgmhOgQP84eiTZEFvuK0DWFEF3iP6jGxR+T0vWTF/0OTBdCVIj/e7EVQhRCu8bxvhCiX/x7bC6EcIj/RqS8BVTiVwxlpxAiFAhBG1oZIKW8nMx+iwALIAgtge154fV+QDRaT/ghWrJ7jpRyPzAI2CGEqJlMjAhgAdqwTRDaeH9XKeVtKWUs0B54D/BDGzrpmcrfsTfaBWF/YBvamPj+ZNp3B+2i7Fdo1xPuAJ+j/X8zQfvW4w/8BzQGRr0i5l/x7XuM9t50kVJGp7K9C4GNwD60v5eVaNcVngKtgF7x7bgPzEEbplPeAkJKtRCLomRFQojvgPeklB9ldluU7EX1+BVFUd4yuiZ+IcR4odUOuSSEWC+EyKVnPEVRFOX1dBvqib+J5ShQRUoZIYTYCPwtpXTSJaCiKIqSKnoP9ZgBFvEzEix5fhqaoiiKkgl0K+wkpbwnhJiPNnMiAtgnpdz34n5CiOFoN+VgZWVVs1KlV93LoiiKoiR15syZICllkbQco+dQTwFgC9pUtCdo9UY2SylTrFleq1Ytefr0aV3aoyiKkh0JIc5IKWul5Rg9h3paAP9KKQPj5x1vRStIpSiKomQiPRO/H1A3vjaLAJqT/J2ciqIoihHplvillCfR6rKfBS7Gx1quVzxFURQldXRdtUdKOZUktVjSIzo6mrt37/Ls2TMDtUpJkCtXLkqVKoW5ufnrd1YUJdvI8su13b17lzx58lC2bFmSrKehZJCUkkePHnH37l3KlSuX2c1RFMWIsnzJhmfPnlGoUCGV9A1MCEGhQoXUNylFeQtl+cQPqKSvE/W+Ksrb6Y1I/IqiKIrhqMT/Gtu2bcPOzu65h4mJCbt3787spr3kyZMnLFmyJLOboShKFpetEn/x4iDEy4/iGVhJtHPnznh7eyc+Ro0aRcOGDWndunWqjpdSEhcXl/4GpIFK/Iry9nC+6EzZRWWhBMktRvRK2SrxP3iQtu1pdePGDaZNm8batWsxMdHeunnz5uHg4ICtrS1Tp2ozV318fKhcuTKjRo3C3t6eO3fusH79eqpVq4aNjQ2TJk1K9vxeXl7Ur1+f6tWrU7t2bZ4+fcqzZ88YNGgQ1apVo0aNGhw6dAiAy5cvU7t2bezs7LC1teXmzZt8+eWX3Lp1Czs7Oz7//HMCAgJo1KgRdnZ22NjY4OHhYZg3QlGUTOV80ZnhO4fjG+ybruOz/HTOpMaNA2/v9B3bpEny2+3sYNGi1x8fHR1Nnz59mD9/PtbW1gDs27ePmzdvcurUKaSUdOjQgSNHjmBtbc3169dZvXo1S5Yswd/fn0mTJnHmzBkKFChAq1at2L59O506dUo8f1RUFD179sTFxQUHBwdCQkKwsLDgp59+AuDixYtcu3aNVq1acePGDZYuXcqnn35K3759iYqKIjY2ltmzZ3Pp0iW849+kBQsW0Lp1a6ZMmUJsbCzh4eHpe/MURclSphyYQnh0+v8/v1GJPzN98803VK1alV69eiVu27dvH/v27aNGjRoAhIaGcvPmTaytrSlTpgx169YFtJ58kyZNKFJEK6DXt29fjhw58lziv379OiVKlMDBwQGAvHnzAnD06FHGjBkDQKVKlShTpgw3btygXr16zJw5k7t379KlSxcqVKjwUpsdHBwYPHgw0dHRdOrUCTs7Ox3eGUVRjM0v2C9Dx79Rif91PfNXzU50d09/XHd3d7Zs2cLZs2ef2y6lZPLkyXz88cfPbffx8cHKyuq5/V5HSpns9MqUju3Tpw916tRh165dtG7dmt9//53y5cs/t0+jRo04cuQIu3btol+/fnz++ef079//tW1RFCXrio6NxiqHFaFRoek+R7Ya49fD48ePGTRoEH/88Qd58uR57rXWrVuzatUqQkO1v4B79+7x8OHDl85Rp04dDh8+TFBQELGxsaxfv57GjRs/t0+lSpXw9/fHy8sLgKdPnxITE0OjRo1wdnYGtGsMfn5+VKxYkdu3b1O+fHnGjh1Lhw4duHDhAnny5OHp06eJ5/T19aVo0aIMGzaMIUOGvPTBpSjKmyUoPIhW61oRGhWKmUn6++1vVI//dYoVS/5CbrFi6T/n0qVLefjwISNHjnxu++TJk+nZsydXr16lXr16AOTOnZt169Zhamr63L4lSpRg1qxZNG3aFCklH3zwAR07dnxunxw5cuDi4sKYMWOIiIjAwsICNzc3Ro0axYgRI6hWrRpmZmY4OTmRM2dOXFxcWLduHebm5hQvXpxvv/2WggUL4ujoiI2NDW3btsXGxoZ58+Zhbm5O7ty5+eOPP9L/RiiKkqkuPrhIhw0dCHgawNrOaxFCMOXAFHxJ+wVe3RZiSY/kFmK5evUqlStXzqQWZX/q/VWUrG/7te18tPUj8ubMy/Ze26n9Tu3E17LaQiyKoihKBkgpmXFkBp1dOlO1aFVODz/9XNJPr2w11KMoipJdhEWFMeivQWy6sol+tv1Y3n45ucxyGeTcKvEriqJkMX7BfnTc0JELDy4wr+U8Pqv3mUGLKqrEryiKkoUc9TtKF5cuRMZG4trblbYV2ho8hhrjVxRFySJ+P/s7zdY0I3+u/JwcelKXpA86Jn4hREUhhHeSR4gQYpxe8RRFUd5UMXExjN09lmE7h9GsXDNODj1JpcKVdIun52Lr16WUdlJKO6AmEA5s0yueXsaPH8+iJLcMt27dmqFDhyY+/+yzz1i4cOFLx9WvXx/Q7uL9888/9W+ooihvpEfhj2izrg2/nPqFCXUn4NrHlQIWBXSNaayhnubALSll+krJpUFCqVKT700ou6gszhedM3S++vXr4+npCUBcXBxBQUFcvnw58XVPT08cHR0Tn8fGxiZuh/Ql/oRzKIqSvV1+eJk6v9fBw8+D1R1Xs6D1ggzdkZtaxkr8vYD1egdJWqpUIvEN9mX4zuEZSv6Ojo6JSfzy5cvY2NiQJ08eHj9+TGRkJFevXiU4OJimTZvSp08fqlWrBmh38QJ8+eWXeHh4YGdnx48//khsbCyff/55YinnZcuWAVo9oBfPoShK9rXz+k7qrqxLWHQY7gPcGWg30Gixdf9oEULkADoAk1N4fTgwHEgsd5yScXvG4X0/5brMJ+6eIDI28rlt4dHhDPlrCCvOrEj2GLvidixqk3L1t5IlS2JmZoafnx+enp7Uq1ePe/fucfz4cfLly4etrS05cuTg1KlTXLp0iXLlyj13/OzZs5k/fz6urq4ALF++nHz58uHl5UVkZCSOjo60atUKIMVzKIqSfUgpmX10NlMOTsG+hD3be22nVN5SRm2DMaZztgXOSimTXQ5FSrkcWA5ayYaMBHox6b9ue2ol9Po9PT2ZMGEC9+7dw9PTk3z58iWO5deuXTtVCXvfvn1cuHCBzZs3AxAcHMzNmzfJkSNHqs+hKMqbKTw6nCE7hrDh0gZ62/RmZYeVWJhbGL0dxkj8vTHQMM+reuYAZReVTXZFmjL5yuA+0D3dcRPG+S9evIiNjQ2lS5dmwYIF5M2bl8GDBwM8V4b5VaSU/PLLLy8t3eju7p7qcyiK8ua5G3KXThs6cTbgLLOaz2KS4ySD3pSVFrqO8QshLIGWwFY94ySY2XwmluaWz22zNLdkZvOZGTqvo6Mjrq6uFCxYEFNTUwoWLMiTJ084fvx4YmXOlLxYKrl169b89ttvREdHA1qp5bCwsAy1T1GUrO34nePUWl6LG49usKP3Dr5s8GWmJX3QuccvpQwHCukZI6m+1foC2rJkfsF+WOezZmbzmYnb06tatWoEBQXRp0+f57aFhoZSuHDhVx5ra2uLmZkZ1atXZ+DAgXz66af4+Phgb2+PlJIiRYqwffv2DLVPUZSsa/W51YzYNYLSeUtzcMBBqhSpktlNUmWZ33bq/VUUfcTExfDF/i/48cSPNC/XnI3dN1LQoqDB46SnLLOq1aMoimJgjyMe02tLL/bd2sfY2mONNj8/tbJOSxRFUbKBq4FX6bihIz5PfFjRfgVD7Ye+/iAjU4lfURTFQP6++Te9t/Qml1kuDg44SAPrBpndpGSp6pyKoigZJKVk7rG5tPuzHeULlMdrmJeuSb94cRBCe0DNmmk9XvX4FUVR0sj5onPi7MFSeUtROm9pPO960r1Kd1Z3XI1VDn3vyXmQ7O2wqacSv6IoShok1AQLjw4H4E7IHe6E3KFb5W64dHMx6vz88tzidjqOU0M9qWBqaoqdnR3Vq1fH3t4+sWiboihvnykHpiQm/aS8/L2MmvQ78BdnSPMoD5BdE38A0Bi4b5jTWVhY4O3tzfnz55k1axaTJydbb85gVFlmRcm6/IL90rTd0KIjYviByfxFJ/7hvXSdI3sm/unAUWCa4U8dEhJCgQLaIgnu7u60a9cu8bVPPvkEJycnDhw4QOfOnRO379+/ny5dugBakbZ69ephb29P9+7dCQ0NBaBs2bJMmzaNBg0asGnTJsM3XFGUDIuOjU5x/N4636urCxuC/7kHXCjeisnMZikf04Cj6TrPmzXGPw5IuSozeABxSZ7/Fv8wARqmcIwd8Orab0RERGBnZ8ezZ88ICAjg4MGDr9y/WbNmjB49msDAQIoUKcLq1asZNGgQQUFBzJgxAzc3N6ysrJgzZw4LFy7k22+/BSBXrlwcPZq+v0hFUfT1IPQB3Td1JzQqFDMTM2LiYhJfM0RNsNfxWuhB6c97UjnuCSMs1rAson+6z5W9evy1gaL8/7cyiX9eJ2OnTRjquXbtGnv27KF///68qtSFEIJ+/fqxbt26xGJubdu25cSJE1y5cgVHR0fs7OxYs2YNvr7/rybas2fPjDVUURRdeN3zotaKWnj5e+HcxRmnTk6UyVcGgaBMvjIsb788wzXBUhIXK9nfej41PmtKpJkVD7afYGl4f6QELQ2dOZPWc75ZPf7X9MwBGIlW3T8XEAV0BZYYrgn16tUjKCiIwMBAzMzMiIv7/1eMZ8+eJf48aNAg2rdvT65cuejevTtmZmZIKWnZsiXr1ydfpVqVZVaUrGeN9xo+dv2YYrmLcWzwMexL2APoluiTenTrCdfqD6Llw+2csu5K1eMrsSqZL8PnzV49foAHwAjgRPyfBrrAm+DatWvExsZSqFAhypQpw5UrV4iMjCQ4OJgDBw4k7leyZElKlizJjBkzGDhwIAB169bl2LFj/PPPPwCEh4dz48YNwzZQURSDiI6NZuzusQz8ayD1S9fn9LDTiUnfGM47nSO0Yk1qP3TFs/uPOPy7ySBJH960Hn9qJK38/6thTpkwxg/aHXpr1qzB1NSU0qVL06NHD2xtbalQoQI1atR47ri+ffsSGBhIlSpaGdYiRYrg5ORE7969iYzUVgWbMWMG77//vmEaqiiKQTwMe0iPTT047HuYcXXGMa/VPKMVWZNxkgN9VtLA5RMemxbm9qrD1B9U36Axsl/i18GrplfOnTuXuXPnJvva0aNHGTZs2HPbmjVrhpeX10v7+vj4ZKiNiqIYxhn/M3R26UxgeCB/dPqDftX7GS12yP1wztYdRQvfNZwr0pJyx5wpUaGIweOoxK+TmjVrYmVlxYIFCzK7KYqipNLa82sZ7jqcIpZFODroKDVLpu8GqfS4tuMGdO9Go6hLHGsxlfq7v0GYmeoSSyV+nZxJ+4V2RVEySUxcDJ/v+5xFJxfRuExjNnbfSFGrokaLf2j0JmouGUKMyMHl+btx/Kz16w/KAJX4FUV5qwWGBdJzc08O+RxibO2xzG81H3NTc6PEjgiO4qjjF7S8/BOX89al6KGNVLMvrXtcXRO/ECI/8DtgA0hgsJTyuJ4xFUVRUutcwDk6uXTiQegDnDo6McBugNFi+3jcIbhND1qGn8DT4VPqHJ6LqUUOo8TWezrnT8AeKWUloDpwVed4iqIoqfLnxT9xXOVInIzj6OCjRk36R7/ZS57GNXg3/BJnv9xI/VOLjJb0QccevxAiL9AIGAggpYxCu6VKURQl08TExTBp/yQWnlhIQ+uGbOq+iWK5ixkldvSzWA42nU7LE9O4bVEVC9fN2DeraJTYSenZ4y8PBAKrhRDnhBC/CyFeujVVCDFcCHFaCHE6MDBQx+akX0JZZhsbG7p37054+MslWTPi559/pnLlyvTt25fIyEhatGiBnZ0dLi4uBo2jKG+7oPAg2qxrw8ITC/nE4RMO9D9gtKQfcCGQcyXa0vrE93hV7Efpeyd5JxOSPuib+M0Ae+A3KWUNIAz48sWdpJTLpZS1pJS1ihQx/HxVQ0io1XPp0iVy5MjB0qVLDXr+JUuW8Pfff+Ps7My5c+eIjo7G29tb1e5RFAPyvu+NwwoHPPw8WNVhFb988IvRLuJ6/eSJrFED2ydH8Bq+gjpXnchZwNIosZOjZ+K/C9yVUp6Mf74Z7YPgjdawYcPEkgudOnWiZs2aVK1aleXLlwOwcuVKxo8fn7j/ihUrmDBhAgALFy7ExsYGGxsbFi3SCg+NGDGC27dv06FDB+bMmcNHH32Et7c3dnZ23Lp1y8i/naJkT+svrqf+yvpEx0bjMciDQTUGGSVuXKxkb5sfsRvXmFjTnARsOY7DsqEJi+VmGt3G+KWU94UQd4QQFaWU14HmwJWMnHPcOPB+VVnmdLCzg0WpKf4GxMTEsHv3btq0aQPAqlWrKFiwIBERETg4ONC1a1d69eqFra0tc+fOxdzcnNWrV7Ns2TLOnDnD6tWrOXnyJFJK6tSpQ+PGjVm6dCl79uzh0KFDFC5cmDp16jB//nxcXV0N+4sqylsoJi6GyW6TmX98Pg2sG7Cp+yaK5y5ulNiPbgdzpf4QWj/YwplSHal0wgmrd/IbJfbr6D2rZwzgLIS4gFb5/ged4+kioVZPrVq1sLa2ZsiQIYA2Nl+9enXq1q3LnTt3uHnzJlZWVjRr1gxXV1euXbtGdHQ01apV4+jRo3Tu3BkrKyty585Nly5d8PDwyOTfTFGyr0fhj2jr3Jb5x+czstZIDvQ/YLSkf2HdBUIq1qLeg+2c6DoPe99tWSbpg87z+KWU3kAtQ50vtT1zQ0sY40/K3d0dNzc3jh8/jqWlJU2aNEksyzx06FB++OEHKlWqxKBB2lfKV9XvVxTFsM7fP09nl87ce3qP39v/zhD7IUaJqxVY+x1Hl7GEmBbg1u+HqDskpVWgMk/2K8tsJMHBwRQoUABLS0uuXbvGiRMnEl+rU6cOd+7c4c8//6R3794ANGrUiO3btxMeHk5YWBjbtm2jYcOs9w9CUd50Gy9vpP6q+kTGRnJ44GGjJf0Q/1A8yvWjhctwrhduQM7L56iYBZM+qJIN6damTRuWLl2Kra0tFStWpG7dus+93qNHD7y9vRPX57W3t2fgwIHUrl0b0L4VvFjGWVGU9IuNi+WrA18x13Mu9UvXZ3P3zZTIU8IosW9suYhZ7+44Rt/kaKtpOO76SrcCawYhpcwyj5o1a8oXXbly5aVtb4IPP/xQurm5ZXYzXutNfX8VJalH4Y9kq7WtJN8hR+wcISNjInWLVaxYwqKHUkKcHMRKGYaFDKCYPL/ooG5xUwKclmnMtWqox8CePHnC+++/j4WFBc2bN8/s5ihKtnfxwUUcVjhw6N9DLG+3nN/a/UYOU/3KHzx4oP1pSRhODGQVQzhOPezwxvbTprrFTVYAVKZymu8CU0M9BpY/f361nKKi6MT5ojNTDkzBL9gP63zWdKzYkZXnVpI3Z14ODzxMvdL1jNKOKlxmE92pxDW+YyrT+YY4MmFoZzpYYpk7rYepxK8oyhvB+aIzw3cOJzxaK5niG+zLz6d+5r0C73F40GFK5ilplHYMwIkljOIpeWjFPg7Qwihxn2MBPEv/4WqoR1GUN8KUA1MSk35SUbFRRkn6YYHhHK0wCCcGcZI62OGdOUkfwAt4J/2Hq8SvKMobwS/YL9ntd0Lu6B775o6r+JeuTf1/1jCNb2jJfu5jnBlDL9kPtAQCAAGStN8kpBK/oihvhAIWBZLdbp3PWreYUoL7kLWU7FiLAlEPOT97D0uKTSM2mVHyYnoX+YwEJgKtgIJAE2AkXOd6mtc5UYk/FRLKMlevXh17e3s8PT0BiIuLY+zYsdjY2FCtWjUcHBz4999/M7m1ipK9RERHMHTHUP6L+A8T8XzKsjS3ZGbzmbrEDXkQgft7Q2myqj+38tdCnvOmxqRW3L///8mcSR/37+vSDM01oC6wABiFNtRzAPgVwgiLSOvp1MXdVEhasmHv3r1MnjyZw4cP4+Ligr+/PxcuXMDExIS7d+9iZfXSkgMZEhsbi6lpFr4RRFF0dPvxbbpt7Ma5++eY0nAKFQtV5JtD3yTO6pnZfCZ9q/U1eNzLW69j2rs7TaMucqzJV9Tb+z0mOTIhXUpgOTAesAT+Ajpk/LSqx59GISEhiXfjBgQEUKJECUxMtLexVKlSia8l5eXlRf369alevTq1a9fm6dOnODk58cknnyTu065dO9zd3QHInTs33377LXXq1OGHH36gR48eifu5u7vTvn17APbt20e9evWwt7ene/fuhIaG6vVrK4rRud5wpebymvz75F929t7JjGYz6Fe9Hz7jfIibGofPOB+DJ30pYe+APynTtSZFY/y5NH83jodmZk7SfwR0AUYADYCLGCTpw5vW48+kuswJ1TmfPXtGQEAABw8eBLSyDA0aNMDDw4PmzZvz0UcfvVSGISoqip49e+Li4oKDgwMhISFYWFi8Ml5YWBg2NjZMmzaNmJgYypcvT1hYGFZWVri4uNCzZ0+CgoKYMWMGbm5uWFlZMWfOHBYuXMi3336bsfdDUTJZbFws37l/xwyPGdgVt2NLjy2UL1Be97iP/SM45TiO1j7LuVLAkRKHN2BTrZTucZN1AOiPtobhAmAcBu2mqx5/KiQM9Vy7do09e/bQv39/pJSUKlWK69evM2vWLExMTGjevDkHDhx47tjr169TokQJHBwcAMibNy9mZq/+vDU1NaVr164AmJmZ0aZNG3bu3ElMTAy7du2iY8eOnDhxgitXruDo6IidnR1r1qzB19dXnzdAUYwkMCyQNs5tmOExg8F2g/Ec7GmUpO/tcp2AsvVo7bMcr2aTqHz/EAUyI+lHAZPQZu3kBU4CEzB4pn6zevyZVZc5iXr16hEUFERgYCBFixYlZ86ctG3blrZt21KsWDG2b9/+XKkGKSUimdV2zMzMiIuLS3yeUNIZIFeuXM+N6/fs2ZNff/2VggUL4uDgQJ48eZBS0rJlS9avX6/Tb6ooxnXy7km6bepGYFig0Uopx8XBnj5/0MhlFFEmubi+wBWHCR/qHjdZN4A+wBnXL3lhAAAgAElEQVTgY2Ah2ri+DlSPP42uXbtGbGwshQoV4uzZs/j7+wPaDJ8LFy5QpkyZ5/avVKkS/v7+eHl5AfD06VNiYmIoW7Ys3t7exMXFcefOHU6dOpVizCZNmnD27FlWrFiRuA5v3bp1OXbsWOIykOHh4apUhPJGklKyxGsJDVc3xNzEHM8hnkZJ+oH/hnLQegAfuAzAr3BNTC94UzEzkr4EVgI1gH+BrcBSdEv68Kb1+DNJwhg/aP9I16xZg6mpKQ8fPmTYsGFERkYCULt27ecu2ALkyJEDFxcXxowZQ0REBBYWFri5ueHo6Ei5cuWoVq0aNjY22NunvByxqakp7dq1w8nJiTVr1gBQpEgRnJyc6N27d2L8GTNm8P777+vxFiiKLsKiwvjY9WOcLzrzYYUPWdt5bYrz9Q3p9Mrz5B/Rk2YxNzj9wbfU3P4NwjwT0uF/wHBgC9AM+IMM3ZGbWkJmoZWhatWqJU+fPv3ctqtXr1K5cuVMalH2p95fJbPceHSDLi5duBJ4helNpzO54eSX5ukbWmyMZG/n32jmOoEQ04KE/ObMe8OMXFEzgTvQD7gPzES7OSsdv74Q4oyUMk0rHer6ESeE8AGeArFATFobpyhK9rTlyhYG/TWInGY52fvRXlq+21L3mPevPeFGoyF8ELiV8yXb8N7RNRQtV1T3uC+JBqYCs4H3gOMYcIHa1DHGGH9TKaWdSvqKokTHRjNx30S6bepG5SKVOTv8rFGS/smfThBV1Y56gTvw6jEPW79dWGVG0v8HcARmAYOBsxg96cMbMsaf0swYJWOy0jCfkv0FPA2g5+aeePh5MNphNAtaLSCnWU5dY8ZExbG/zXxaHJrCA/NS3HU+ikOvOrrGTJYE1gBj0LLuJqBbxk+7bl36jtO7xy+BfUKIM0KI4cntIIQYLoQ4LYQ4HRgY+NLruXLl4tGjRypJGZiUkkePHpErV67MboryFjjiewT75facCTjDus7rWPzBYt2T/r1zDzlT7APaHprEhXKdKOhzjnKZkfSfAL2AQUBN4AIZTvoxMTBxIvTrl77jdb24K4QoKaX0F0IURSsmOkZKeSSl/ZO7uBsdHc3du3efm+euGEauXLkoVaoU5ubmmd0UJZuSUrLw+EImuU3i3YLvsqXHFmyK2uge13PGQcp/25f88jEXBy/C4fePITNGDTyAjwB/4Hu0m7MyWHrr8WPo1Qv27YPRo+HXX7PYxV0ppX/8nw+FENuA2kCKiT855ubmlCtXTo/mKYqio5DIEAb/NZgtV7fQpXIXVndcTd6ceXWNGRUeg3uz72lxcia+OSsSuWkvDu1tdY2ZrGhgGvADUA44hpb9MujyZejYEfz8YMUKGDoUfv017efRLfELIawAEynl0/ifW6G9FYqiZHOXHl6i68au3PrvFvNbzmdCvQm6X6fzPXaXR2360CrUgxOVB1Lj6GJyFjRstdxUuQ30BU4AA4GfgTwZP+1ff8FHH4GVFbi7Q/366T+XnmP8xYCjQojzwClgl5Ryj47xFEXJAv68+Cd1fq9DSGQIBwcc5LP6n+me9D0m7SRPw+q8H3qW05+upe6V1ZmT9NcBdsBVYAOwmgwn/bg4mD4dOnWCSpXg9OmMJX3QsccvpbwNVNfr/IqiZC1RsVF8tvczFnstpqF1Q1y6uVAij77LE0Y8icSz0Zc0v7iIG5Z2WO50oVazTLh7PRhtgZQ/0UoorwPKvPKIVAkNhYEDYcsWrbe/fDm8prhvqqhaPYqiZNid4Ds0Wt2IxV6LmVhvIgf6HzB40i9eXLs+m/CoLK5yvUAdml9cxDH7MZS7f5xSxkz6AUBjYCdaL98FbTD7EAZJ+rdvaz37bdtgwQL44w/DJH14Q+bxK4qSdbnddqP3lt5ExkSyuftmulbpqkucBw8SfpIM5Xd+4lPCsKIdO3E9006XmK/0PdqsnSNoF3A9gHqGOfXBg9C9uzbMs3s3tGplmPMmUD1+RVHSJU7GMfPITFqtbUUxq2J4DfPSLeknyM9jNtGdFQzHk/pU5zy7MHLStwAEsAztTiXQqmo2y/ippYSff9YSffHi4OVl+KQPqsevKEoqOV90ZsqBKfgF+/FO3ncobFEY7wfe9KnWh+XtlmOVQ9+LqQ3wwJm+lCCAL5jDfCYijd13jQO+RBvSSVhOwxLoDMzP2KkjI2HkSFi9WpuyuXYt5DHAbKDkqMSvKMprOV90ZvjO4YRHhwNwN+Qud0PuMqD6AFZ3XK3rrJ3oiBg8Wk7DnZn8Sznq48lpHHSLl6J/0OrreADWwB0gJ/AMbbWs4uk/dUAAdOkCJ07At9/C1KlgouNnmkr8iqK81pQDUxKTflLuPu66Jv3bB3142rEvzUI9WUN/PmExoYaYFJ8WccAStLtuzdGmaO4A2qHV0l+OdqE3nU6e1JJ+cDBs3gxd9R0tA1TiVxQlFfyC/dK0PaOkhIMfu1BzxccUIY5TnzozaUMfQh+8vG+xYro0QfMvWi/fHWgDrABKod2YlSAdd84mWLMGhg+HkiXB0xNsjXSTsbq4qyjKK53xP4OpSfIFZqzzWRs8XpBPKAfLDqb5il4E5K1E2DFvai/qw/372gfCi4/79w3eBK2X/xtQDW0N3N+Bv9GSvgHExMD48doc/QYNtIu4xkr6oBK/oigpkFKy6MQi6q2sR94ceclp+nw1TUtzS2Y2n2nQmJ6LzxL8nj1N/Zw40XwKFR96ULx+eYPGeC1ftAIzo4D6wCVgCNpMHgN49AjatIFFi+DTT2HvXihc2DDnTi2V+BVFeUlQeBAdNnRg/N7xfFDhA26OvcnKjispk68MAkGZfGVY3n45fav1NUi8Z+FxbG+0gFpj6mIlwrm94gB13WZgktOIlWMl2lCODXASbbrmXrQLuQZy6RLUrg0eHrBqlZb8zTJhwD3Lr7mrKIpxHfE9Qp8tfQgMD2R+y/l8UvsTXS/gXj10n/86DMAxdB/ny3fi/cO/Y1GqkG7xknUHGArsQ5uPvxIoa9gQ27Zp9fPz5oWtW6FuXcOcNz1r7qa6xy+EKCOEaBH/s4UQwsiX1hVF0VNsXCzTDk+j6ZqmWJpbcnzIccbUGaNb0o+Lg7+G7qRQM1vsQ49w+ZPfqP7PVuMmfQmsQuvlH0O7ULsfgyb9uDj47jtt5k7VqlqRNUMl/fRK1ZcMIcQwtIlLBYF30S5xLAWa69c0RVGMxf+pP3239sXdx52PbD9iyQdLyJNTv75dwD9hnGnyGR3vLeN23uqYuTpTtWFV3eIl6x4wDNiNVnNnFWDgywlPn0L//rB9OwwYAEuXQlZY9C61Pf7RaEsEhwBIKW8CmbBSsaIohrb75m6qL63OqXuncOroxNrOa3VN+ofmnSasoj0f3FvO+ZYTKffgJAWNmfQT1r+tijZN82fgIAZP+rduQb16sHMn/PijdkduVkj6kPp5/JFSyqiEr3xCCDP+X6VCUZQ3UFRsFF8d+IoFxxdgW8wWl24uVCpcSbd4ocGx7G8+m3ZnvuOReXHurXKj+gADFLhJC3/gY8AVrXzyauA9w4dxc4MePbSf9+yBFi0MHyMjUpv4DwshvgIshBAt0SY67dSvWYqi6On249v02twLL38vRtUaxYLWC8hlpl931Hv7v8T07kfnZ8e4ULknlQ79Ro5iBXSL9xKJVit/DBAB/Bj/cwbXv30pjNRm6kycCJUra6tmvfuuYWMYQmoT/5doM1kvon1e/o12S4OiKG+YjZc3MmznMEyECVt6aOvh6iUmWuLaax3Nto4GIbg2ZS220/sad+Hz+8AI4C+0sslOgA5l+589g48/1urmd+6s3ZWrV5G1jEpt4rcAVkkpVwAIIUzjt71cvENRlCwpPDqccXvGseLsCuqWqsv6ruspm7+sbvF8vR/zT4sRdHq0kWtFGlDSbS2VbPWL9xKJtjjKaCAMmAeMx+C9fIB797RZO6dOwfffw9df61tkLaNS27QDaIk+gQXglpoDhRCmQohzQgjXtDZOURTDuPzwMrVX1GbF2RV86fglRwYe0S3pSwn7Jx/E1N6WRo+24t19JpUC3MlrzKT/EOgO9EYbw/cGJqJL0j9+HGrVgitXtLn6336btZM+pL7Hn0tKGZrwREoZKoSwTOWxn6ItPZw3rY1TFCVjpJSsPLeSsbvHkidnHvZ+tJdW7+qwske8x/cjOdrkaz68voC7FhUI3HAcuw5purco4zahXYUMAWYDn6FbOcrVq2HECChVCvbvBxsbfeIYWmo/l8KEEPYJT4QQNdEukbySEKIU8CHqeoCiGF1IZAi9t/Rm2M5hOFo7cn7EeV2T/qnVl/G3rkP76/Pxrv0x79w/S0ljJv0goBfQA+0GrLNopZR1SPrR0TB2LAweDI0aaUXW3pSkD6l/S8YBm4QQ/vHPSwA9U3HcIuALSLmAthBiONrNYVhbG77Sn6K8jU77n6bn5p74PvHlh2Y/MKnBJEyEPuMPUc/i2NNuMa0OfEGYaV5uLNiJ/QQjL4e4De0C7mNgJlrW0amXHxSkTdU8dEirsDl3bubU28kQKWWqHmhLENigFSo1T8X+7YAl8T83AVxfd0zNmjWloijpFxsXKxd4LpDm08xl6YWl5VHfo7rGu3nAV57M00xKkOdLfyhDb93XNZ6UUkp/KWUjKWWAlDJIStlHahmkhpTygr6hz5+XsmxZKXPmlHLNGn1jpRZwWqYyjyc8Xvk5JYRoJqU8KIR4cb5XBSEEUsqtrzjcEegghPgAyAXkFUKsk1J+lPaPJ0VRXicoPIgB2wfw982/6VSpEys7rKSgRUGDnLt4cXjw3CIokn6s5RfGUJw4zo1cTo1fhxpnmuZ04CgwCO2ibRDwPTAZrXuqk82btbIL+fPDkSNalc031eu+oDRGu5m5fTKvSSDFxC+lnIz2V4EQogkwUSV9RdHHYZ/D9Nnah6DwIH5p+wujHUYbtLha0qRfmECW8TFd2MYRGlLpuBM16hqhZr4F2vq2CfbE/5kD+Fa/sHFx2hq4M2ZoxdW2boUSJfSLZwyvTPxSyqlCCBNgt5Ryo5HapChKKsXGxTL9yHSmH5nOewXfY1efXdgVt9MtXgf+YgXDyEcwE5nHj4wntq4OcySTcxvtyqJH/HMzoBvaXbg6CQnRSinv2KFdyF2yBHLmfP1xWd1rL0lIKeOEEJ8A6U78Ukp3tHJIiqIYyL2Qe/TZ2ocjvkfoZ9uPJR8uIXeO3LrEykswixjHIJw4hx3NOMhljDiN5QEwlv8n/RxADFAAKK5PyJs3oWNHuHEDfvkFRo827g3Hekrttej9QoiJaPfBhSVslFL+p0urFEV5pV03djFg+wCexTxjTac19K/eX7dYx384xAUGUoq7TOdrpvMN0eTQLd5zJLAObV5hGFAFaAiMBJYDAfqE3bsXevUCU1PYtw+aGbmWnN5Sm/gHo/0VjHphu5EXw1SUt4/zRWemHJiCX7AfpfOVpkrhKuy5tYfqxarj0s2FioUr6hL3v3sRnGrxFW2uLeIGFXDkGCcx4goid9Aqg+1GW/t2JZC0eOivhg8pJSxcCF98oc3L374dypUzfJzMltrEXwUt6TdA+wDwQFuIRVEUHTlfdGb4zuGER2tlsfyC/fAL9qNluZbs6LNDt4qah+d7UfLL/rSJvcZJh0/o5TsHn4cv36xfrJgOwePQevNfALHAT2j1dnS+lBARAcOHw7p10LUrODlBbn1GzjJdahP/GrQboH+Of947flsPPRqlKIpmyoEpiUk/qRv/3dAl6T+6H417y5l0vDSDIPMS/LN4P3VGtOBfg0dKwU20VbEOAy3QPgCM0OO+exc6dYIzZ2D6dJgyJfuM5ycntYm/opSyepLnh4QQ5/VokKIo/+cX7Jem7Rnh9uNFCn8xiK4xZzhfvR+V9/1M8aL5DR4nWbFos3O+AXKiFXkZDBgh+R47pvXww8O1+vkdOugfM7Ol9h7uc0KIxME9IUQdtKWJFUXRyZYrW1Kci2+dz3DlTQL9o9lgO5NGE2pijR8+8zdT3fsPchgr6V9Cq5P/OdAKuIK2+ocRkv6KFdC0qVY3/8SJtyPpQ+oTfx3AUwjhI4TwAY4DjYUQF4UQF3RrnaK8hR5HPKbftn5029SN0nlLvzSkY2luyczmMw0Sa9+Ci/iXqUuvi19z06YLefyuUPazrgY592tFod1xaw/4ABuA7UBJ/UNHR8Mnn2hj+k2banX0q1TRP25Wkdqhnja6tkJRFAD23drH4L8Gcz/0PlMbT2VKwylsvLIxcVaPdT5rZjafSd9qfTMU56F/DAdazaHr5e8JNcuP74LNVJ1gpIQP4MX/1/Trg3YBt7BxQgcGQvfucPgwfP45zJqlTdt8q6S1uI+eD1WkTXlbPY18Kke6jpR8h6y8uLL0uuelS5y4OCl3z7soz5nWlBLk5Wo9ZLT/Q11iJStcSvm5lNJESvmOlHKH8UJLKeW5c1KWKaMVWVu71rix9YKhi7QpiqK/o35HGbh9ILcf32ZC3QnMaDYDC3OL1x+YRvfvxnCg9Ry6X9F6+XcWbqLK+G4Gj5OiI2i9/H/QZu7MA/IZL/zGjTBwIBQqBEePaqtmva2y+AJhipJ9PYt5xhf7v6DR6kbEyljcB7qzoPUCgyd9KWHXnEvcL1uXvle+5rZtZ/L6Xaa0sZL+U7S7gBqjzd45gDZN00hJPy5Om57ZsyfUqKEtmvI2J33QbakCRVFe5VzAOfpt68flwMsMtx/O/FbzyZMzxfWK0i3gjjaW3+Oa1su/t2gTlT41Yi9/D9oyS3fRyi7MAKyMFz44GD76CFxdYehQWLw4exRZyyiV+BXFiGLiYpjlMYtpR6ZRxLIIf/f5m7YV2ho8jpTgOvsSpb4ZyEexZ7hevQfv7VlMweJFDB4rWf8B44E/gMqAJxiz2gNoxdU6dIBbt+DXX2HkyOx9U1ZaqMSvKEZyLega/bf1x8vfi942vVn8wWKDLZSS1L3bkRxuM4tuN38gLL6XX9GYvfzNaCUW/gO+jn8YuZe9ezf07g3m5uDmBo0bGzd+VqfG+BVFZ3EyjkUnFlFjWQ1uP77Nxm4b+bPrnwZP+lKC6zcnCa5Qkz43v+e2fXfy3rnMO8ZK+veBrkB3oBTalM3pGDXpS6mtgfvhh1C2LJw+rZJ+clSPX1F05PPEh0F/DcLdx51277djRfsVFM9t+ALy926Ecar113T0+YmgHO/gv9iVSsM+NHicZEm0IZ3xQDgwC5iI0bNLeLg2jr9+vbYY+qpVYGXE6wlvEpX4FUUHUkpWnVvF+L3jAVjZYSWD7AYZdDlELQ7smehGlR+H0Vn6cKHBKGx2zsIkf16DxkmRH1rp5D1oq2yvBPSpEv3qZvhB585w7hz88AN8+aUaz38VlfgVxcACngYwbOcwdt3cRZOyTVjdcTVl85c1eJw7Fx9zpc1ntPVfjZ/F+9xbdQTbXg0NHuc5AUAvYD1aeYVJaD3+n9HG9TNh8NjDA7p108oq79gB7doZvw1vGt0SvxAiF9otGznj42yWUk7VK56iZAUbL29k5K6RhEeHs6j1IsbUGYOJMGw2lBLcRm7BdvlomssgzrSaTI1t32JiqU9t/udMR1uNoxbah0ALYAVQVv/QyVm2TKu5U768VoKhUqXXH6Po2+OPBJpJKUOFEObAUSHEbinlCR1jKkqm+C/iP0b/PZoNlzZQ+53arOm0hkqFM56FiheHBw+SPCeAxXxCV7ZyM3cNYtbvpma7GhmO81oWwLMkzxOWPDxKpiT9qCj49FNYuhTatoU//4T8Riommh3o9sUsvoxEaPxT8/iH1CueomSWv2/+jc0SGzZf2cyMpjM4NviYQZI+JE36ksGs5ApV+IC/mcRs3n10ineMkfQBnIGklw0sgL5gvBVa/u/hQ2jRQkv6kybBzp0q6aeVriNyQghTIYQ38BDYL6U8mcw+w4UQp4UQpwMDA/VsjqIY1NPIpwzfOZwP//yQQpaFODX0FFMaTcHMxLBfpCtyjUM0ZSVDuYAt1TnPXCZhksMIl+juA/3QpmnGodXIz4X2fT4vYPgJSq909qxWbsHLS+vlz579FlbWNABdE7+UMlZKaYc2q7e2EMImmX2WSylrSSlrFSlipLsKFSWDDvscxnapLSvPrWSS4yRODztNjRKG7X0/uf+MqXzHeapTnfMMZQVNOcRN3jdonGTFAovRFjd3AaYAzYCRwAlgBNqHghFt2AANGmg/Hzum3aClpI9RZvVIKZ8IIdzR6vpfMkZMRdFDRHQEUw5OYdGJRZQvUJ4jA4/gaO1o0BhSwv6vDlF+3gi+4wbr6MtnLOAheqxsnoyTaEXVzqJdvF3My1M0fzVOUwBiY7Uia3PmaIl/yxYoWtR48bMjPWf1FAGi45O+Bdo/oTl6xVMUPThfdE5cBKVY7mKYYIJ/qD+jao1iTss55M6R26DxrhwJwrf7RNo+XMOdnO/SMnYfbrQ0aIwU/QdMRpulUxxtRaweGGUJxJQ8eQJ9+mglGEaMgJ9+ghw5Mq892YWePf4SwBohhCnakNJGKaWrjvEUxaCcLzozfOdwwqPDAbgfeh+BYJLjJGa3mG3QWCHBkl091tBy30QqEIz3B19hu/FrLr5rAQ9e3r+YITv/ccAa4AvgMVoVze94/mJuJrh2DTp2hNu34bfftMSvGIbQFnDJGmrVqiVPnz6d2c1QFADKLiqLb7DvS9vL5CuDzzgfg8SQEnYtvE6Br0bgGOXOP8XqU3jLcvI7VjXI+V/rPNqwjifanbdLAFvjhH6VXbu0nn7OnNrQTkOd70t7kwkhzkgp07TCgCrSpijJeBzxONmkD+AX7GeQGFe9I1lT/ntaTrSlWqw3/365jPf8PYyT9EPQauvUBG4Aq9Fut8zkpC+ltgZu+/bw3ntakTWV9A1PJX5FSUJKyR/n/6Di4pQLzljns85QjNBQWNHTDZMatgz0+Q6/Wl2x8r1KuVnDwUTn/5ISbey+EtoC50OB68BAMj0bhIVpM3W++gp69dJKMVhn7K1WUqASv6LEuxp4lWZ/NGPA9gG8W/Bdfmj2A5bmls/tY2luyczmM9N1filh17K7HCrag2EbW5I/TyxPNuyhgtefmL5jhAnx14GWQG+0K3AngKWA4ZcESDNfX23GzsaN2uwdZ2ewtHz9cUr6qCJtylsvPDqcGUdmMN9zPrlz5GZZu2UMtR+KiTDBOr914qwe63zWzGw+k77V+qY5xo1LURzu8hO9b36PmYjFb9h0rH+eCLmMUF8nHJiJtri5Jdr0zBFAFrnx6fBhrchadLQ2tt/W8AuSKS+SUmaZR82aNaWiGJPrdVdZdlFZyXfIAdsGyAehDwx6/rAwKVf0OSivUFlKkLdsOsjoG7cNGuOV/pJSlpHa/7D+Usr7xgv9OnFxUi5ZIqWZmZQVK0p5/Xpmt+jNBJyWacy1qsevvJXuBN/h0z2fsu3aNioXroz7AHcalzXcUk1Swt7V/kSOmcjQ8PU8zF2Ox0t2Ur6fkWoG/wt8CuwEqgDuQBZaiSoqSququWKFtlqWszPky5fZrXp7qMSvvFWiY6P5+eTPTHWfSpyMY1bzWUyoN4Ecpoa7K+jWtWgOdllMz6tTySmi8B04lTJLJoGFhcFipCgSmA/MQBvKmYf2AWCuf+jUevAAunbVyi589RVMm6bq7RibSvzKW8PzjicjXEdw8eFF2r3fjl/a/mLQBVIiImDDqCM4rBnNMHmJ25U/oPTWnylT6V2DxXglN7TFUG6gFVX7EShtnNCpdfq0tlLWo0da7Z2ePTO7RW8nNatHyfYehT9i2I5hOK5y5MmzJ2zruY0dvXYYNOkfcLqDW5HeDHJqTDHLp/y3ajvlL7tiboyk74+2KlZLtOJqu4HNZLmk7+yszck3NQVPT5X0M5Pq8SvZlpQSJ28nvnD7gscRj5lYbyJTm0w1aH0d32sRHOs8j07XZmMiJP/2+5ZySycZZy5iDPALMBWIQiuzMAmtbHIWEhurrYE7fz40bgybNoEqxJu5VOJXsqVLDy8xctdIjvodxbG0I799+BvVilUz2Pkjn0lcB27GwWUiffDjWrXulN8yj3IVyhgsxisdQyu1cAFoi/YBYKQRpbR4/Fi7KWvvXhg9Gn78Ecyz0PWGt5Ua6lGylbCoMCbtn0SNZTW4EniF39v/zpFBRwya9D1/O493gaZ0delBXL4CPNzoTqULG8mhZ9IPQJuVcwkYDDRAK6i2FdhFlkz6V65A7dpw8CAsXw6LF6ukn1WoHr+Sbey4voMxu8fgF+zHYLvBzGk5h8KWhdN9vhfXuy1EEDP4mmGsIMSkAJfHLKXqj0ONMyXle7RFzmuiVdOcBHwDWOkfOj127oS+fbURr0OHwNGwSxYoGaR6/Mobz/eJLx03dKTjho7kyZEHj0EerOy4MkNJH/6f9M2IZiw/cZMKDOV3fmEMlvduUvXnj/VP+hZo9fCXodXZiUIb2/+JLJn0pYSZM7Vyyu+/ry2RqJJ+1qMSv/LGio6NZu6xuVRZUgW3227MbTGXcx+fo4F1AwNFkHyIKxew5SfG4YUDtlxgPIvIWbyAgWKkGBr+AsrFP09YDMWSTFvk/HVCQ6FHD/j6a6237+EBpbPYzCJFo4Z6lDeSh68HI3eN5HLgZTpV6sRPbX7KcNXMBFKC56/ncGMizTnIDSrQke3soAO6L0clgb1owzingffQ1rp1B3ICz8iURc5f599/oVMnuHRJm70zYQKITFy5S3k11eNX3iiBYYEM+msQjZwaERoVyo5eO9jWc5vBkv7p7XfZU2Ig9cbUpDrn+YRfqMpldtAR3ZP+IbSLtm2BQGAVcBXIh1ZULZMWOX+dQ4fAwQH8/ODvv+Gzz1TSz+pUj1/JkpKudWudz5oZzWbwLOYZk9wmERIZwiTHSXzT6BuschhmoPvyiadc6j+H9jcXYkos3i0+p5nbZILJb5Dzv9IxtB7+ITZguZMAAB5uSURBVOAd4De0mTsJVSS2JtnXiIucv46U2kyd8eO18fy//oIKFTK7VUpq6LnYemngD7QvpXHAcinlT3rFU7KPF9e69Q32ZcD/2jvzsCir9o9/DiiIu+CCioAbWu4KaWXupm/uZT+zMjML99J6W8wyK31bbBdTc83MJV8NzSxzt95SwF1UXEEFBcVYRPY5vz/OKCCMsczwAHM+1zXXMDPPPM/Ncfh65j73+d6BIzFJEw95PsS8fvNoUds6XarOncpg14jF9AuazjBiCG09nIar/kP7e72p4A7xtux3G4IS/F+B2sAXwBhK3AasvEhNhfHjYckSGDgQvvsOqhrco1eTf2w5488AXpFSHhBCVAH2CyG2SimP2/CamjLAtO3Tbov+LUzShJuLG7uf3Y2wQh7hymXJev9f6bbp3zzHcc57dCZ+yU+06H1f1jG2SqkcAaajFm9dgY9QHjslsEonLy5fViZrf/0Fb78NM2bYvnGYxrrYTPillJdR206QUiYKIU6gvshq4dfcFUs9ba8nXy+y6MfFwerJe7nnu6mMN+0ipmoTYj9dT8PRg22fmD6BslX4AbVA+x7KObMUzZSDgpTJWlycsl4YOtToiDSFoVhy/EIIb6AdsC+P1/wBfwBP3WDTrknLTGPxgcUIIVD9JXJSlAXcmzdh1dvHqTNnGmPTA4l3rk3Ma3Oo/ZY/OFnPkjlPzqI2YH2PqsufBrwC2Lgi1NosXw7+/lC3rprttza4Mbum8Nhc+IUQlYF1wGQpZcKdr0spvwG+AfD19c39164p82SaMll1bBXv7HqHc3+fw8fVh4j4CFIzU28fU9het+np8MOnF3B8bwbPJn9LqmMlosa/T72PJkNl65m15ckF4H1gKcoP/2XgNaCUGZRlZMDrr8Nnn0H37qovbs2i7Y3TGE1BW3YV5Ib6uG8BXs7P8br1on1hMplk4IlA2fLrlpIZyDbz2sifT/0sTSaTXHFkhfT63EuKGUJ6fe4lVxxZUaBzZ2ZKuW7BVbm42hSZgpNMEc7y4v+9LOXVqzb6bbIRKaWcIKV0Mt8mSSmjbH9ZWxAbK2Xv3lKClJMmSZmWZnREmjuhEK0XbSn6AlXV80V+36OF337YcW6H7LSok2QGsulXTeXqo6tlpimzyOc1maTcsi5Rfu3+royniszAQUb0fk6awiOsEPU/EC2lfFlKWUFKWU5K6S+lvGD7y9qKY8ekbNRISicnKRcvNjoajSUKI/y2TPU8CIwAjgohDpmfe1NKudmG19SUcIIjg5m2Yxpbz23Fo6oHCwcsZGSbkZR3LLpt494dN9n/wnyGnvuIh4nhgu8QKi2dhWfLe6wQ+V24jmp3+BWQjPrUTwca2faytiQwEEaMUNmwXbvg/vuNjkhjTWxZ1fMHNt/qqCktHL96nLd3vs36E+upWbEmnz38GeP8xlGhXNGL1o+FpPDns98wIPQDOnGFi816kr5oFp6dO1oh8rsQj6q9/wxIBIahqnaa2faytsRkgpkz4Z131G7cH3+E+vWNjkpjbfTOXY1NCY8LZ8auGXx35Dsqla/EjK4zmHL/FKo6F72G8dyJVPaMXEzv4Fn4E0VEw64kz1tDgz5drBD5XbgBBAAfozzxH0VV7bS07WVtzY0bMHIkrF8PzzwDCxZAhVKwmUxTcLTwa2xC9I1oZu6ZyYL9C3AQDkzpNIU3Or9RZKtkgCsX0tj5zFIe3D2LZ7nIeY/OJASswGtQdytEngeXUT1tlwGBwAcoL51+qFr89ra5bHFy7pyyUj5+XHXJeukl7bdTltHCr7EqcSlxzP7fbL7Y9wWpGamMbjeat7u+jUdVj3yf484GKLeoVyudOX7L6fDL+wyXEZyr3YlrXyym4RO9bKtSM1BNUFoCN4FeKMEvI3nv7duVnbKU8Ouv0Lu30RFpbI0Wfo1VSEpLYk7QHD7630fEpcTxRMsneK/bezR1K7hr152i70Qqz7CcN65+SOPN5zjr6kfU7Pk0GtXHtoJfAUjN9viWi8QflAnRlxK++kq5aTZvrkzWGpfAFo4a66OFX1Mk0jLTWLh/ITN/n8mVG1fo17QfM3vMpK172yKf24WbvMBCXmU2HkQShB/iiy9o/GJ/2wr+MVSHq1s4oGwGKwJDUBU8pZyUFBg3DpYtUz76y5dDlSpGR6UpLrTwawpFpimT749+zzu73iE8LpyHPB9i7eNrrdL9qgoJTGAuU/ic2lxlF10ZxVK20Qv5ko0E3wRsRlXpbEdZK4wE4oD/omb/JbQJSkGJioJHH4V9+1T1zvTp2mTN3tDCrykQUko2hG3grR1vEXo1lHbu7Zj31Dz6NO5TZAO10D2xnH3xSy7wFdWJ5xf6Motp/A9rtVLMg0TUou1XwBnAA7V4+wLghqrYGYtyk/oGs+1g6WXvXiX6CQmqemfIEKMj0hiBFn5Nvtl+bjtv7niToMggfNx8+GHoDzx272M4iMJPF00m2LY8irjpn/HIxfm0IIn1DGEW0zhABytGfwfngTnAYiAB6ATMRAl99r1kJbQJSmFYtgzGjFF1+Vu2QKtWRkekMQot/Jpc3Nn9alTbUfx+4Xe2n99Og6oNWDRgESPbjqScQ+E/PjduwMYPQnGa8ykDE1fgSCYn2w3HI2Aq4x9tkWdVT5EboEhgDyp/vwGVu38cZY1s471eRpKRAf/+N3z5JfTsCWvWgJub0VFpjETIPOxvjcLX11eGhIQYHYZdc2f3q1tULl+Z93u8z1jfsUXabRsRLtn8+m4arZ9Nn4zNpDi4cLHXc3h/OYXyzW1UUpICrEYJ/iFUCmcMMB7VIaIMExurSjV37IDJk2H2bCinp3tlCiHEfimlb0Heoz8Cmhy8uf3NXKIPUMOlBpM7TS7UOaWEv37PYN/r6+i89xPGEUK8cy0ujnqPBh+Mp6mtpp9XgPmoHrYxQAtgIfAUavG2jHP0qNqUFRkJS5fCs88aHZGmpKCFXwNAcnoyK4+utNj96lLCpQKfMz0dAlfcIGLGEh678DlTCOdqDR+uv7YA15dGUM3FRup7ADW7XwWkA/1R6Zye2I171Lp1yn6halXYswc6luFUlqbgaOG3c87/fZ55IfNYfHAx15OvU96hPOmm9FzHFaT7VWwsrPnwPA7z5zLsxmJqEMflRg+SMutzav3fQNvUDmai8vZfoHbZVkJV40wCCr6HrNRiMsG778J77ymxX78e6tUzOipNSUMLvx1ikia2ndtGQFAAm05twkE4MOSeIUz0m8ilhEv4b8qZ489v96sTxyW/vr6TJpu/YqxpI1I4EP3QY1T7z0vU7fyAbX6ZOFRlzhwgAvAGPgWeA6rb5pIllcREZaW8YQOMGgVff61N1jR5o4XfjohPiefbw98yN3gup2JPUbtSbaY9NI0xvmNyeukIclT1zOo5i6daPZXnOaWEbRtvcmLaCnqEfsUUQkl0duPaM1OpPX0c9Tzy79FTIE6hau+XAUlAV+BzYCDgaJtLlmTOnFH5/LAwZcMwcaI2WdNYRlf12AGhMaHMDZ7L8sPLSUpPopNHJyb6TWTovUNxLudcqHPevAmBX4ST/OnXDLm+CFf+5op7WypOfZGqLzwB1sjf33LFXIPaLSuBbaj8/c+AEzAclb9vV/TLlVZ++w2GDVMZtLVroUcPoyPSFCe6qkdzmwxTBhvDNhIQFMDO8J04OzozvNVwJvhNwLdegT4jOYiMyGDHq7/gHjifJ9J/wYQDkX5DqPLhi7h372zdaeb7KEO0t4H7UIIfCtRGOWaOBYpa21+KkVJZKL/6KrRoobpmNSrFXb80xYcW/jJGTFIMiw4sYl7IPC4lXMKzmicf9vyQ0e1H58sL35IlcovqkbzbYBH3HV3ECC7xt7M7l0a8SYP3/fHyyv/Cb75wQdXe32KR+SaAb1Gdrgr3RaXMkJysduF+9x089pjalVu5stFRaUoLNhN+IcQSVCFdjJSylPcmKvkERQYREBTAmtA1pGWm0atRLwL+FUB/n/44OuQ/6Z1d9B3I5GF+YwwL6B+3iXJxmZz0epgrr3+J+/MDqFG+6H1yc3EJeA21WPu3+TlHoAdK9Ota/5KljchI5bETHKyqd6ZN0yZrmoJhyxn/MlSDuuU2vIZdk5KRwg+hPxAQFEBwVDCVnSrj396f8X7juadW4RuMe3CRZ1jOCyzEmwiiqc1sXmXioRdo3sYGuYTrKAfMVcBuVC6/FmqG74SqxW+CFn3gzz/VDP/GDZXaGTTI6Ig0pRFbNlvfI4TwttX57ZmL8ReZHzKfhQcWcvXmVZrXbE7AvwIY0WZEoXvZXjyVzOF3A9nCUnqxDQck2+jJq8xmA4NIx4mpbaz4SyQBPwErgV9R4t4MlbsfDryOEvoy4oppDRYvVh76np6wbZvK62s0hUHn+EsJUkp2he8iIDiAwJOBAAxsNpCJfhPp0bBHoSyRoyIlv3+yD6dVy+gevZr+xBOOF+8xneU8w3msPLtPB35Dif0GlPjXR1XlPAm0JWtnbRlyxSwq6enw8ssQEKDaIq5eDa6uRkelKc0YLvxCCH/UvA5PTysvEpZC7nTGnN51OqkZqQQEB3D86nHcXNx49YFXGec7Dq/qXgU+f3Q0/LIoktRF39E1fBnDCCPZoSJn2g0l6eVnaTSiKxIrJoxNqMqcVcBaIBZwBZ5GiX1nsOblyhrXrsHjj8OuXapF4ocfapM1TdGxaR2/OdWzKb+Lu/Zex2/JGROgfd32TLpvEsNaDMOlfMFq5K9dg5+/u861BetoH7aSruzGAckFz86UHzOKupMev913z1JVT506cOVKPi8ogcOomf0q1IJtRWAwKo3zMCp3r7krhw+rHP6VK7BoETz9tNERaUoiuo6/lPPGtjfyFH33Su6EvBBSoHTO33/DpjVJRM7bSIsjqxjOrziRzjVXH6498Q61Jz+JZ9PcJjb5Fve8OIMS+pXASdSnqy/wMWpHbaUinNvOWLtWuWnWqAG//w5+fkZHpClL2LKccxXQDagphLgEvCOlXGyr65VWzlw/w4aTGwgMC7TogBmdFJ0v0U9IgE3r0zj99VZ8QlbyqAykEjeJr1KfuMEvUuulJ6nZvp11N1ldBn5AiX2Q+bmuwBTgMZT3vSbfmEyqB+6sWfDAA8pl072U9/jVlDxsWdUz3FbnLs1IKTlw+QCBJwMJDAvkWMwxANrUaUM152rEp8bnes/dnDGTkuCX9cmcnvsbnsHreMT0E08SR1IFV5L+NYKKLw6nWpeHrFvoHYdafF0J7ETl8dsDs1GbqxpY71L2REKCSuf89BOMHg1z54KznW9U09gGneopBtIz09kdsZvAk4FsCNvApYRLOAgHunh14fM+nzOo2SAa1miYZ44/L2fM5GTYuj6RcwGb8QheT9/MnxlKEknONUjqMQjT2KFU6vswlZyKkEi/0ycnGdiEEvvNQBqqtv4tVN6+eeEvpYHTp1U+/9QpVb0zfrw2WdPYDi38NiIxNZEtZ7cQeDKQn0//TFxKHC7lXOjTpA8zu8+kn0+/XBYKtxww83LGTE2FHWtjCZ/7M57B63g4cwsVSCXBpTYJvZ7GZfxjVOrZjUrW2k17yyfneVQVzo/ADVRt/QRURU4H7KaxiS359VcYPlxV62zbBt26GR2Rpqyj3TmtSPSNaH469ROBJwPZdm4bqZmpuLm4MaDZAAY3G0zvxr2pWL6ixffnrqiR3MtxhpTbxCNyEx0z/8QRE9creXCj96PUe/ExynV5EByt6ENcAUjN43lHVA1+V+zS9tgWSAmffAJvvAGtWqmduN7eRkelKW3oqh4DOB17+na+/q+LfyGReFf3ZpzvOAY3H8yDng9SziF/wxwdDc6k0JXd9GcT/dlEQ8IhAyJc23G+zzS8J/bH9X4/XK2VB0hE2SRsN9/uFH0n4FGU171eZLQaycnw/POwcqWq01+6FCrpqidNMaGFv4CYpIn9Uftvi/3xq8cBaOfejhndZjC4+WBa1W6V79LL+DjJwbVniP/vVn7kN3qxjcokcRMXttGLD5jKz/QjMra+dX6BVGAvWUK/D9W2sAJqM9WTQDAQiBL9NKAGWvStyMWLymTtwAFVvTN1qs7na4oXLfz5IC0zjV3hu9hwcgMbwjYQmRiJo3Cki1cXxnQYw6Bmg/K9izY6GvZuvs71tdupGrSVDrG/0Y0IAMLxYjnPsIn+7KQ7KVihmUkmcIgsof8dtVDrAPihPHF6Ag+gxB/UDH8s2ifHBvzxhzJZS06GjRuhf3+jI9LYI1r4zdxplfBWl7eo6lz19uJsQmoCFctXpG+TvgxqNoh+TfvhVvHuRepSQng4/Lk9maj1e6mybysdrm9lAPtxQJLkWJULzXsQ1uc1PJ/rTcM2TSjyaqlEtSW8JfQ7ybI3bgG8gBL6rkA1C+fQPjk2YeFCmDBB5fF37YJ7Cm+gqtEUCb24y92tEmpWrMlAn4EMbj6YXo163dUuwWSCEyfgr603uLrhTyrv3027xN3cRxBOpJMpHLns1Ql69cZ9RG/KPXBfDuOVu33dv+s/UxRK5LcBO1AWCQCeKJHvifKz17bGhpCWBpMnw7x50LcvrFoF1e2sEbzGdujF3QIQlRhFSFQIwZHBfPLXJ6RkpOQ6pk6lOkS+HGmxkUlGBhw8CPu2xBG36Q+qHd7NfSl7eJb9lCOTTOHIVW9f4rtNxm1IFxy7dcGjqmXb5Dp1LPvk5OBvYBdZs/qT5ufdUAJ/S+wbo8stDSYmRi3e7tkDr70G//mPdYuwNJrCYBfCH3szVol8VDDBUcGERIUQlRgFgINwwCRNeb4vJikmh+gnJ0PQnxmErQ8leedeXE/vxTdjLxPNypvu4ESsT0cSe71B9UFdcXzgftwL0A/vtk9OXpuntpIl9AdQu2UrAl2A0UAvoDXa6bIEcfAgDB6sxP/77+HJJ42OSKNRlDnhT0hN4MDlAwRHZon8+bjzt1/3cfOhu3d3fOv54lfPj7buban6ZgtMVSJynUskeLJ9xWWiftwH+/bRIGovvjKYriQBkFihJomtO5HY+2mq9H2Q8h074u5ihQXZd1GLsLeMzf5EVdeUAzqhmo/3BDqiXS5LKGvWwKhR4OamFnQ7dDA6Io0mi1It/DfTb3LoyqHbs/mQqBDCroUhUQlx7+re+NbzZazvWHzr+dKhbgeqVci9omn6bRb0f4GGScm0uwLtL0O7KAfaR8ThnlEPgAxRjst12nHV9zkY2InKPTtRpWFDqhS1Di8OOGq+vYiqwrlFsPneEWWT8BCgG2qXaDIz4a23lG9+587w3//mkarTaAymRC3uinpCer3mddumIDtpmWkcjT6aI2UTGhNKplRK6V7ZHb96fupW348OdTtQq1KtXNeIjYXTR5KJ2XOSlJBjVDh5iCpnDtDWIYgaJrW4myEg1NmTgyndeXBCOxo86keF+9tBUWbz6UAYSuCPZLu/mO2YKqiSyr+BDMAFVVr5CbqOvhQQH6/SOZs3g78/zJkDRbFL0mjyQ5lY3I2Ij8D/J38uxV+iTuU6t1M2h6MPk5aZBoCriyt+9fwY4DMAv3p++NbzpX7VrA1Oyclw5gz8EZpK7P9OknE4FJdzodSKCaVpeih+nMMRlddPERU4TGtWm57hIO04QHuOyZakpqiidhlQwF9Aoqps7hT4EyjxBzXqzVEbploDrcz3HsB4VO38LeuEqmjRLwWEhSmTtbNnVfXO2LFGR6TRWKbEzfgZk/O5Kk5V6FCvA751ffGrr2b03tW9kVJw4QKcOmkiMiiSGwdPYwo7jUvkadwTTuHDKZpwhnLm3EmmcCSmWlMSvVrCvS2ocn8LanVtQbl7fRDlLf//d9fhuQGEklPgjwLXsx1Tn5zi3gol+pZmgo+Su8n4egvHakoEmzcrkzVnZ5Xa6dLF6Ig09kRhZvwlV/glnJh4Ajd8OHssjai9F4g7HEFqWDjOF8/gev00jU2nacIZXMgqxUxzcCauZhNSvXwo36YF1R9sQYUOLcDHx6K5+a00vTuwGmUpf6uqUkpU3v0suQX+bLaTVEKJenaBb4VyttSUSaSEjz9Wlgtt2iiTNa+Ct0HWaIpEqRf+am5CPt0JvOLA62pFGp9rRf3MCOqSsx9guijP9RqNSfFoimPzplRp35SqHZoifJqCh0eBm47ccsWcC4xB2c7vBDpWgOEtUbP6ZPPBDkBTcgp8a8AbXUppR9y8qZqlrF4Nw4bBkiVQ0bLxqkZjM0q98PsKIUOAFEe44OhOZq2WZHh44dTEi2qtvajp60W5Rl7QoEH+dsFkotIu17LdrubxeAtgaRh6klPg7wVrWOhoSi8XLqj6/EOH4IMP1MYsbbKmMYoysLjbjKXNJzC6XHXksRHI7C1oJcpC+BqwH8sinv3xdSwLemWgpvnWFQhHWR1kAM5AP9RXAL2wqsnGnj0wdCikpsKmTfDII0ZHpNEUHJsKvxCiL/AlqhJ9kZTyw7u/ozKjTk5iFOZy9p7kFPI0C28rhxLwWub71tl+rnnHazVR1gZ3ztrHkVVNkwbUQYu+Jgfz58OkSdC4MWzYAM2aGR2RRlM4bCb8QghH1Jy5N2ouHSyE2CilPP5P7001v6FxCip37ktu8c7+uCpF96SJRlsRW4l/yh7mJ7tYko5JT4dXXoEFC9QMf+VKqGbJ2VSjKQXYcsZ/H3BGSnkOQAixGhgEWBT+g6gMTDqQLqDi4ZyvF8sf+pJsx1SwfFhJEqb8HGOt69gzU6fC++9rkzVN6cdmi7tCiKFAXynl8+bHI4COUsqJdxznj5pn44hjqwY0iC9P+fKnOX0210nti5qoBJe9o8chCz0WWeixyKKZlLJKQd5gyxl/XsmXXP/LSCm/QSVXEEKEnJfnC7Q6XVYRQoQUdKW+LKLHIQs9FlnoschCCFHgJia2rDy/BDTI9tgDZWag0Wg0GgOxpfAHA02FEA2FEE4oh/mNNryeRqPRaPKBzVI9UsoMIcRE1PYoR2CJlDL0H972ja3iKYXosVDocchCj0UWeiyyKPBYlKiduxqNRqOxPdpdRqPRaOwMLfwajUZjZxgi/EKIJUKIGCHEsWzPuQohtgohTpvvaxgRW3FjYSxmCyFOCiGOCCF+FEJUNzLG4iKvscj22r+FEFIIUdOI2IobS2MhhJgkhAgTQoQKIT42Kr7ixMLfSFshxF4hxCEhRIgQ4j4jYywuhBANhBA7hRAnzJ+Bl8zPF0g/jZrxLwP63vHcG8B2KWVTYLv5sT2wjNxjsRVoKaVsDZwCphZ3UAaxjNxjgRCiAcr640JxB2Qgy7hjLIQQ3VG731tLKVugmnLaA8vI/bn4GHhXStkWmG5+bA9kAK9IKe8BOgEThBD3UkD9NET4pZR7yNmnCtQH+lvzz98Cg4s1KIPIayyklL9JKTPMD/ei9kCUeSx8LgA+B17DstdqmcPCWIwDPpRSppqPiSn2wAzAwlhIlEsXQDXsZI+QlPKylPKA+edEVFPX+hRQP0tSjr+OlPIyqF8OqG1wPCWF54BfjA7CKIQQA4FIKeXhfzy47OMDPCSE2CeE2C2E8DM6IAOZDMwWQlxEffOxl2/FtxFCeAPtgH0UUD9LkvBr7kAIMQ311e57o2MxAiFERWAa6qu8Ru27qYH6iv8q8IMQdtsCZhwwRUrZAJgCLDY4nmJFCFEZWAdMllImFPT9JUn4o4UQdQHM93bxNdYSQoiRQH/gKWm/my0aAw2Bw0KIcFTK64AQwl47JVwC1ktFEGBCmZXZIyOB9eaf16LcgO0CIUR5lOh/L6W8NQYF0s+SJPwbUf+YmO83GBiLoZgb2LwODJRS3jQ6HqOQUh6VUtaWUnpLKb1RwtdeSnnlH95aVgkEegAIIXwAJ+zXoTIK1TsP1JicNjCWYsP8DW8xcEJK+Vm2lwqmn1LKYr8Bq1CtTtJRf8yjUX2xtqP+AbcDrkbEVkLG4gxwEThkvs03Ok6jxuKO18OBmkbHaeDnwglYARwDDgA9jI7TwLHojGrCehiV4+5gdJzFNBadUQvbR7LpwyMF1U9t2aDRaDR2RklK9Wg0Go2mGNDCr9FoNHaGFn6NRqOxM7TwazQajZ2hhV+j0WjsDC38Gs1dEEJMM7sgHjE7QXY0OiaNpqjYrPWiRlPaEULcj9o93V5KmWq2hHYyOCyNpsho4ddoLFMXuCaz3DDtdZespoyhN3BpNBYwG2H9AVQEtgFrpJS7jY1Koyk6Osev0VhASnkD6AD4A1eBNUKIZw0NSqOxAnrGr9HkEyHEUGCklHKA0bFoNEVBz/g1GgsIIZoJIZpme6otEGFUPBqNtdCLuxqNZSoDc8zN7jNQrqn+xoak0RQdnerRaDQaO0OnejQajcbO0MKv0Wg0doYWfo1Go7EztPBrNBqNnaGFX6PRaOwMLfwajUZjZ2jh12g0Gjvj/wGgxCyJPAY+tAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(S,price_0, color='blue', marker='s',label=\"Zero costs\")\n",
    "plt.plot(S,price_w, color='green', marker='o',label=\"Writer\")\n",
    "plt.plot(S,price_b, color='magenta', marker='*',label=\"Buyer\")\n",
    "BS.plot(axis=[10,20,0,8])  #plot of the Black Scholes line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec5.1'></a>\n",
    "### Time complexity\n",
    "\n",
    "If we set the \"Time\" argument to True, the method also returns the execution time.\n",
    "Let us verify that the algorithm has time complexity of order $\\mathcal{O}(N^4)$\n",
    "\n",
    "The following operation will be very time consuming. In case you are in a hurry, reduce the NUM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>N</th>\n",
       "      <th>Price</th>\n",
       "      <th>Time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>50.0</td>\n",
       "      <td>2.241215</td>\n",
       "      <td>0.033321</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>100.0</td>\n",
       "      <td>2.249143</td>\n",
       "      <td>0.096112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>200.0</td>\n",
       "      <td>2.245423</td>\n",
       "      <td>0.395083</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>400.0</td>\n",
       "      <td>2.246784</td>\n",
       "      <td>1.948342</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>800.0</td>\n",
       "      <td>2.246289</td>\n",
       "      <td>13.928330</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1600.0</td>\n",
       "      <td>2.246577</td>\n",
       "      <td>111.881915</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>3200.0</td>\n",
       "      <td>2.246412</td>\n",
       "      <td>1607.674912</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        N     Price         Time\n",
       "0    50.0  2.241215     0.033321\n",
       "1   100.0  2.249143     0.096112\n",
       "2   200.0  2.245423     0.395083\n",
       "3   400.0  2.246784     1.948342\n",
       "4   800.0  2.246289    13.928330\n",
       "5  1600.0  2.246577   111.881915\n",
       "6  3200.0  2.246412  1607.674912"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "NUM = 7\n",
    "price_table = pd.DataFrame(columns=['N', \"Price\", \"Time\"])\n",
    "for j,n in enumerate([50 * 2**i for i in range(NUM) ]):\n",
    "    price_table.loc[j] = [n] + list(TC.price(n,Time=True))\n",
    "display(price_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the computational times we can estimate the exponent $\\alpha$ of the polinomial growth $\\mathcal{O}(N^\\alpha)$. \n",
    "\n",
    "For higher values of N, the exponent converges to the expected value of $\\alpha = 4$.\n",
    "\n",
    "Here we are quite close."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The exponent is:  3.8449269495846137\n"
     ]
    }
   ],
   "source": [
    "print(\"The exponent is: \", np.log2(price_table[\"Time\"][6]/price_table[\"Time\"][5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec5.2'></a>\n",
    "### Is the risk aversion important?\n",
    "\n",
    "The coefficient $\\gamma$ measure the risk aversion of the investor. We can see how the option price is affected by this coefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt4lPWZ//H3TRIIxwBJUA6BAAIJBBRBC+Ju66n2sGptrbVSrYeW2sqvRd291tXd9tdWd9v9FWu3trLYrtYu29paa61trafSrRuwgkvlEBAV0CBKEjkfQiD3749n8mQymSSD5JlJZj6v63ouZub5zsz9kGTu+Z7N3REREQHok+kARESk51BSEBGRkJKCiIiElBRERCSkpCAiIiElBRERCSkpiIhISElBRERCSgoiIhLKz3QAx6ukpMTLy8szHYaISK+yevXqencv7apcr0sK5eXlrFq1KtNhiIj0Kma2LZVyaj4SEZGQkoKIiISUFEREJNTr+hSSaWpqora2lsOHD2c6lIwqLCxkzJgxFBQUZDoUEemlsiIp1NbWMnjwYMrLyzGzTIeTEe5OQ0MDtbW1jB8/PtPhiEgvlRXNR4cPH6a4uDhnEwKAmVFcXJzztSWRbLRs7TLK7y6nz1f7UH53OcvWLovsvbKipgDkdEJoof8DkeyzbO0yFvx6AQebDgKwbc82Fvx6AQDzp8/v9vfLipqCiEg2ajrWxN89+XdhQmhxsOkgtz9zeyTvmTU1hUy76aabGDduHIsWLQLgwgsvpKysjB/84AcA3HLLLYwePZqbb765zfPOOussqqur2bp1K9XV1Vx55ZVpj11EMmvP4T1srN/IxvqN1NTXhLdf3fUqR5uPJn3O63tejySWnEwKy9Yu4/Znbuf1Pa8ztmgsd5535wlXw8466yx+/vOfs2jRIpqbm6mvr2fv3r3h+erqau6+++7w/rFjx8jLy6O6uhqArVu38l//9V/HlRRaXkNEej535429b4Qf+PFJ4K39b4XlCvoUMKl4ElUjqrhs6mUsWbWEhkMN7V5vbNHYSOLMuaQQVfvcvHnzuOmmmwBYv349VVVV7Nixg127djFgwABqamrYs2cP55xzDiNHjmTNmjVs2LCBQYMGsX//fm699VZqamo47bTT+PSnP80Xv/hFbr31VpYvX05jYyM33ngjn/vc51i+fDlf/epX27yGiPQcjUcb2fzO5uADv66GjQ3Bh/+m+k0caDoQlhtaOJTKkko+eMoHqSipoLKkkoqSCsYPG09+n9aP5srSyjafWQADCgZw53l3RhJ/1iWFRU8sYs1bazo8v7J2JY3HGts8drDpINf/6nruW31f0uecdvJp3P2Bu5OeazFq1Cjy8/N5/fXXqa6uZu7cuWzfvp0VK1ZQVFTEjBkz6Nu3L3/+859Zt25du2Gj3/jGN/jWt77F448/DsDSpUspKirihRdeoLGxkXnz5vH+978foMPXEJH0eefQO60f/PUb2dgQ3N6yewvN3hyWG1c0joqSCv7q9L8KP/grSioYMXBESoNDWr6sdnfrRkeyLil0JTEhdPX48Zg3bx7V1dVUV1dz8803s337dqqrqykqKuKss84C4Mwzz0zpw/zJJ5/kpZde4uGHHwZgz549bN68mb59+6b8GiJyYpq9mW27tyVt7687WBeW65fXjyklU5g1ahbzp88PvvmXVjJp+CQG9h14wnHMnz4/siSQKOuSQlff6MvvLmfbnvaLBY4rGsfya5af0Hu3dBqvXbuWqqoqysrKWLx4MUOGDOG6664DYODA1H5B3J3vfve7XHjhhW0eX758ecqvISKpOdR0iJcbXm7zob+xfiObGjZx+Gjr3J+SASVUlFRwyZRLqCxt/dY/rmgceX2yo38v65JCV+48787I2ufmzZvH4sWLmTBhAnl5eQwfPpzdu3ezfv167rvvPtatW9fhcwcPHsy+ffvC+xdeeCH33nsv5557LgUFBbz88suMHj36hGMUyVXuTt3BunadvBvrN7Jt9zYcB8Awxg8bT2VJJedPOD/84K8oqaBkQEmGryJ6OZcUomyfmz59OvX19W1GEE2fPp39+/dTUtL5L9OMGTPIz8/n1FNP5ZprruFLX/oSW7du5fTTT8fdKS0t5dFHHz3hGEWy3dHmo2zdvbVde//G+o28c+idsFz//P5UlFQwd8xcrj3t2rC9f1LxJArzCzN4BZll7p7pGI7L7NmzPXGTnZqaGiorKzMUUc+i/wvJFfuP7GdT/abWb/6xjt7N72zmyLEjYbmTBp4UftOP7+gtKyqjj+XO/F0zW+3us7sql3M1BRHpPdydt/a/lbSj9429b4Tl8iyPicMnUlFSwYcnfbhNk8+w/sMyeAW9j5KCiGRc07EmXtv1WruO3o31G9nTuCcsN6jvICpKKnhv+XvbfOufOGwi/fL7ZfAKsoeSgoikzZ7De9jUsKndN/9X3nmlzXIOowaPorKkkk/N+FSbpp9Rg0dp4ceIKSmISLdyd7bv2560o/fNfW+G5fL75DNp+CQqSyq5tOLS8IN/SskUhvQbksEryG1KCiLyrhw5doTNDZvbdfRuatjE/iP7w3JD+g2hsqSSCyZc0KbJZ8KwCRTkaZfAnkZJQUQ6tevQrqQdva/teo1jfiwsN7ZoLBUlFVxXdl1rk09pJScNPElNPr2IkkI3ycvLY/r06bg7eXl53HPPPeHSFiI9XbM388aeN9p19NbU17DzwM6wXN+8vkwunsypJ5/KFVVXhB/+k4snM6jvoAxegXSXyJKCmZUBDwInA83AUnf/TkKZ9wG/ArbEHnrE3b8WVUxt7ACuAB6KRXiC+vfvz5o1wUJ8v//97/mHf/gH/vjHP574C3dAy2YLHP8y8IeaDoUreMZ/899Uv4lDRw+F5YYVDqOytJKLJl/UpqO3fGh51iznIMlFWVM4Ctzi7i+a2WBgtZk95e6Jaz3/yd3/JsI4kvs68BzwNeD73fvSe/fuZdiwYGz08uXL26x+unDhQmbPnk1ZWRn33HMPv/zlLwF46qmnuPfee3nkkUd48skn+cpXvkJjYyMTJ07k/vvvZ9CgQZSXl3Pdddfx5JNPsnDhQq644oruDVx6lc6Wgb9w4oVJO3q37NrSZjmH8qHlVJRUcE75OW0md5UMKFGTT46KLCm4+w6C7+O4+z4zqwFGA9FuALAI6HjlbPgTQb2lxb2xow/wVx085zSg83X2OHToEKeddhqHDx9mx44dPPvss52WP/fcc7nxxhupq6ujtLSU+++/n2uvvZb6+nruuOMOnn76aQYOHMg3v/lN7rrrLr785S8DUFhYyHPPPdd5MJITbnvmtqTbNF79y6vbLN1cmF/IlOIpnDHqDK6ecXX4zX9S8SQGFAxId9jSw6WlT8HMyoGZwPNJTs81s78AbwJ/6+7rkzx/AbAAYOzYE9xt6EzgNaCeIDn0AUqAiSf2svHNRytWrODqq6/udAE8M+Oqq67iP//zP7n22mtZsWIFDz74IE888QQbNmxg3rx5ABw5coS5c+eGz/vEJz5xYoFKr+Pu7Dywk3U717Fu5zrW161n3c51HW7H2OzN3PX+u8KO3rFFY3NqOQc5MZEnBTMbBPwCWOTuexNOvwiMc/f9ZvYh4FFgUuJruPtSYCkEax91+oZdfKMH4POxVysEjgAfo1ubkObOnUt9fT11dXXk5+fT3Nz6re3w4dZleK+99louuugiCgsL+fjHP05+fj7uzgUXXMBPfvKTpK+tZbOz265Du8IP/fgEUH+wPixT3L+YqhFVDOo7qM3QzxbjisZx09yb0hm2ZJFIk4KZFRAkhGXu/kji+fgk4e6/NbPvm1mJu9cnlu1WbwM3ENQ9lhJr5Oo+Gzdu5NixYxQXFzNu3Dg2bNhAY2Mjhw8f5plnnuHss88Ggt3aRo0axR133MFTTz0FwJw5c7jxxht55ZVXOOWUUzh48CC1tbVMnjy5e4OUjDpw5AAb6jaEH/7r6taxfud6tu/bHpYZ3HcwVSOq+MiUj1A1oio8WnbsSuxTgGi3aZTcEOXoIwN+CNS4+10dlDkZeNvd3czOJGjMab9DdXeLT0/f656XbOlTgKC6/6Mf/Yi8vDzKysq4/PLLmTFjBpMmTWLmzJltnjd//nzq6uqYOnUqAKWlpTzwwAN88pOfpLEx2A3ujjvuUFLopRqPNrKpYVPrh3/s2LJ7S1imML+QqaVTOW/CeUwrnRZ++JcNKeu0szfd2zRKbohs6WwzO5ugW3ctrV27twFjAdx9iZktJGjMOQocAm529+rOXjfbls5euHAhM2fO5Prrr++W1+vN/xe92dHmo7z6zqttvvmv27mOzQ2bwwle+X3ymVI8haoRVW0+/CcMm6BhnhK5jC+d7e7PAZ2OaXP3e4B7ooqhp5s1axYDBw5k8eLFmQ5FUtSyZ298u/+6neuoqa8J1/A3jInDJ1I1oorLKi9j2oggAUwunkzfvL4ZvgKRzmlGcwatXr060yFIB9ydHft3BJ29O9e3afc/0HQgLFc2pIyqEVW8f+L7wxpAZWmlhnpKr5U1ScHdc36yTW/bRa+naDjY0G60z7qd69h1eFdYZsTAEVSNqOL6mdeHzT5TS6dSVFiUwchFul9WJIXCwkIaGhooLi7O2cTg7jQ0NFBYmLt7y3Zlb+PeNiN+WhLAW/vfCssU9SuiakQVl0+7PPzwn1Y6jdKBpRmMXCR9siIpjBkzhtraWurq6jIdSkYVFhYyZsyYTIeRcYeaDrGxfmO74Z7b9mwLywwoGMC00ml84JQPUFXaOtxTm7hIrsuKpFBQUMD48eMzHYakWdOxJja/s7ndcM9Xd70aLvNQ0KeAytJK5o2dx4LSBeGHf/nQcs3yFUkiK5KCZLdjzcfYsntLmw7fdTvXsal+E03NTQD0sT5MGj6JGSfN4MrpV4ZDPk8Zfoo2chE5DkoK0mO4O7V7a9sN99xQt6HNss7lQ8upGlHF30z6m6DNf8Q0KkoqKMxXf4rIiVJSkIxoWeAt8dv/3sbW5bFGDhpJ1Ygqbph9Q9jhO7V0KoP7Dc5g5CLZTUlBIrX78O7wgz++BlB3sHVQwPD+w6kaUcWnpn+qdcTPiGkM7z88g5GL5CYlBekWB44coKa+pt1wz9q9tWGZQX0HMa10GhdPubjNAm/aw1ek51BSkONy5NgRNtVvajfc87Vdr4U7evXL68fU0qm8r/x9bYZ7lhWVacSPSA+npCBJHWs+xqu7Xm033HPzO5s52nwUgDzLY0rJFGaNmsXVp17dZoG3/D761RLpjfSXm6VS3dC92Zt5fc/r7Tp8a+pqaDwWLN1tGBOGTaBqRBUfrfxoONxzcvFk+uX3S/eliUiElBSyUEcbuu8+tJvJxZPb7ewVv3vXmCFjqBpRxfnjzw87fCtLKhnYVzu+ieSCyPZTiEqy/RSkrXF3j+tw/94WpQNK23T2tizwNrRwaJqiFJF0yvh+CpI+9Qfreb72eVbWrmTl9pWdJoRnr36WaSOmMWLgiDRGKCK9hZJCL3Pk2BFeevulIAHEjld3vQoEHb8zTprR6Ybu54w/J90hi0gvoqTQg7Us+xAmgO0reXHHixw+ehgIZvzOGTOHBbMWMGfMHGaNnMXAvgO1obuIvGtKCj3IgSMHWL1jdZgEnt/+PG/uexMIxv7PGjWLL8z+AnPGzGHOmDmMGTIm6aQvbeguIu+WOpoj1tHQ0GZvZnPD5vDDf2XtSl56+6Vwk/eJwyaGH/5zxsxhxkkztL+viLxrqXY0KylEKFkzTkGfAipLKnlj7xvhdo+D+w7mPWPew5zRQQI4c/SZ2ulLRLqVRh/1ALc9fVubhADQ1NxETX0N15x2De8Z/R7mjJlDRUkFeX3yMhSliEgrJYUINBxsYMmqJby+N/nQ0KPNR1l60dI0RyUi0rXIkoKZlQEPAicDzcBSd/9OB2XPAFYCn3D3h6OKKWqvvvMq3175be5fcz8Hmw5SmF8YjhSKN7ZobAaiExHpWpRLVh4FbnH3SmAOcKOZTU0sZGZ5wDeB30cYS6Sq36jmYz/7GJO+O4mlq5dy+bTLeemGl/jBxT9gQMGANmU1NFREerLIagruvgPYEbu9z8xqgNHAhoSi/wf4BXBGVLFE4VjzMR7d+CiLVyxmRe0KhhUO49azb2XhmQsZNXgUANNPmg5oaKiI9B5p6VMws3JgJvB8wuOjgUuBc+klSeHAkQPcv+Z+vr3y27y26zXGDx3Pdz/4Xa457RoG9R3Urvz86fOVBESk14g8KZjZIIKawCJ335tw+m7g7939WGc7b5nZAmABwNixmWmP37FvB/f8+R7uXXUvuw7vYs6YOfzr+f/KRyo+opFDIpI1Ip2nYGYFwOPA7939riTntwAt2aAEOAgscPdHO3rNdM9TWLdzHXetuItla5fRdKyJSysv5Za5t3BW2Vlpi0FE5ERlfJ6CBV/9fwjUJEsIAO4+Pq78A8DjnSWEdHF3ntnyDItXLOaJV55gQMEAPnv6Z1k0ZxGnDD8l0+GJiEQmyuajecBVwFozWxN77DZgLIC7L4nwvVOSuATFV8/5Kn3ow+IVi/nL23/hpIEnccc5d3DD7BsoHlCc6XBFRCKXs8tcJFuCwjAcZ2rpVG6ZewtXTr+SwvzCE34vEZFMy3jzUU93+zO3t1uCwnFGDBjBus+vS7r6qIhItoty8lqP1tHuZHUH65QQRCRn5WxS6GipCS1BISK5LGeTwvUzr2/3mJagEJFcl5NJofFoIw+tf4ih/YZSNqQMwxhXNI6lFy3V7GMRyWk52dH8z3/6Z9bXrefxTz7Ohyd/ONPhiIj0GDmTFOLnJDjOvLJ5SggiIglyovnoC7/5Alc9chXb9mzDCeZlvLjjRZatXZbhyEREepasTwrL1i5jyaolYTJocejoIW5/5vYMRSUi0jNlfVK4/Znb2yWEFh3NVRARyVVZnxQ6++DXnAQRkbayPil09MFvmOYkiIgkyPqkcOd5d7bbJ9kwbph9g+YkiIgkyPqkMH/6fJZetJRxRePCSWo//uiP+f6Hv5/p0EREepycXTpbRCSXpLp0dtbXFEREJHVKCiIiElJSEBGRkJKCiIiElBRERCSkpCAiIiElBRERCSkpiIhISElBRERCkSUFMyszsz+YWY2ZrTezLyUpc4mZvWRma8xslZmdHVU8IiLStSi34zwK3OLuL5rZYGC1mT3l7hviyjwDPObubmYzgJ8BFRHGJCIinYispuDuO9z9xdjtfUANMDqhzH5vXXxpIHSwG46IiKRFWvoUzKwcmAk8n+TcpWa2EfgNcF064hERkeQiTwpmNgj4BbDI3fcmnnf3X7p7BfAR4OsdvMaCWJ/Dqrq6umgDFhHJYZEmBTMrIEgIy9z9kc7Kuvt/AxPNrCTJuaXuPtvdZ5eWlkYUrYiIRDn6yIAfAjXuflcHZU6JlcPMTgf6Ag1RxSQiIp2LcvTRPOAqYK2ZrYk9dhswFsDdlwAfA642sybgEPAJ7227/oiIZJGUk4KZjQMmufvTZtYfyI+NKkrK3Z8DrLPXdPdvAt9MNQYREYlWSs1HZvZZ4GHg32MPjQEejSooERHJjFT7FG4kaA7aC+Dum4ERUQUlIiKZkWpSaHT3Iy13zCwfTTQTEck6qSaFP5rZbUB/M7sA+Dnw6+jCEhGRTEg1KdwK1AFrgc8BvwX+MaqgREQkM1IdfdQf+A93vw/AzPJijx2MKjAREUm/VGsKzxAkgRb9gae7PxwREcmkVJNCobvvb7kTuz0gmpBERCRTUk0KB2LLUABgZrMIZiCLiEgWSbVPYRHwczN7M3Z/JPCJaEISEZFMSSkpuPsLZlYBTCFYumKjuzdFGpmIiKRdp0nBzM5192fN7KMJpyaZGV0thy0iIr1LVzWF9wLPAhclOeeAkoKISBbpNCm4+1fMrA/wO3f/WZpiEhGRDOly9JG7NwML0xCLiIhkWKpDUp8ys781szIzG95yRBqZiIikXapDUq8j6EP4QsLjE7o3HBERyaRUk8JUgoRwNkFy+BOwJKqgREQkM1JNCj8i2GDn32L3Pxl77PIoghIRkcxINSlMcfdT4+7/wcz+EkVAIiKSOal2NP+vmc1puWNm7wH+J5qQREQkU1KtKbwHuNrMXo/dHwvUmNlawN19RiTRiYhIWqWaFD4QaRQiItIjpLog3raoAxERkcxLtU/huMUmuv3BzGrMbL2ZfSlJmflm9lLsqDazU5O9loiIpEeqzUfvxlHgFnd/0cwGA6vN7Cl33xBXZgvwXnffZWYfBJYS9F+IiEgGRJYU3H0HsCN2e5+Z1QCjgQ1xZarjnrISGBNVPCIi0rXImo/imVk5MBN4vpNi1wO/6+D5C8xslZmtqqur6/4ARUQESENSMLNBwC+ARe6+t4My5xAkhb9Pdt7dl7r7bHefXVpaGl2wIiI5Lso+BcysgCAhLOtolzYzmwH8APiguzdEGY+IiHQuytFHBvwQqHH3uzooM5Zg97ar3P3lqGIREZHURFlTmAdcBaw1szWxx24jmA2Nuy8BvgwUA98PcghH3X12hDGJiEgnohx99BxgXZT5DPCZqGIQEZHjk5bRRyIi0jsoKYiISEhJQUREQkoKIiISUlIQEZGQkoKIiISUFEREJKSkICIiISUFEREJKSmIiEhISUFEREJKCiIiElJSEBGRUO4khR3Ae4G3Mh2IiEjPlTtJ4evAn4DTUWIQEelApNtx9gj9gcNx93cAI4FC4FBGIhIR6bGyv6bwGsmv8jBBwhARkVD2J4WRwPyEx/Jij21JfzgiIj1Z9icFgP3ANFo3Bz1G0Hx0csYiEhHpkXIjKTwCTAY+D/xH7LGnMheOiEhPlf0dzS0eibu9hWA00k+BKzITjohIT5QbNYVEXwbmAp8DtmY2FBGRniQ3k0I+sCx2++PAX6O5CyIi5GpSABgPLAFWAc8BX8tsOCIiPUFkScHMyszsD2ZWY2brzexLScpUmNkKM2s0s7+NKpak+gNXxm47cC/B6CTNXRCRHBZlTeEocIu7VwJzgBvNbGpCmXeALwLfijCO5F4jSAqJSWAm8DjQnPaIREQyLrKk4O473P3F2O19QA0wOqHMTnd/AWiKKo4OjQSGAI0Ecxb6EHQ+vwFcRDCv4T7aLpEhIpLl0tKnYGblBN/Bn3+Xz19gZqvMbFVdXV33BfY2cAOwMvbvyQQ1iGXAAGABMBb4KtCNbysi0lOZu0f7BmaDgD8Cd7r7Ix2U+b/Afnfvshlp9uzZvmrVqu4NMhkHlgOLgd8Q1CY+DdwETIn+7UVEupOZrXb32V2Vi7SmYGYFwC+AZR0lhB7LgHMI+hc2AJ8CHgAqgEuA/yZIHCIiWSTK0UcG/BCocfe7onqftKgk6F/YRjDx7X8INuw5k2BW9NHMhSYi0p2irCnMA64CzjWzNbHjQ2Z2g5ndAGBmJ5tZLXAz8I9mVmtmQyKM6cScRNC/8DrBENY9wCeBU4BvA/syF5qISHeIvE+hu6WtTyEVzQTNS98i2NVtCMHSGV8ExmQwLhGRBD2iTyHr9QEuJuhf+DPwQeAugtnSnwL+N1ZO+0OLSC+hpNBdziDoX3gFWAj8imA/6POAz6ClNESkV1BS6G7lBP0LbxAsvPcs8FuCpqaWpTT6Eizf3bta7kQkBygpRGUoQYf0FQRJAIL/bSOYvz0BKAU+APwT8BjwZvrDFBGJlzub7GTCSILkcJRg8tsR4LMEM6VfIFih9QXgXwi2CAUYRdAUdQYwO3YUpzVqEclhSgpRa1lKYwGwlKDTeVbsaHGQoFO6JUm8QNAn0WICrUniDIK+isFRBy4iuUhJIWrx87i/10GZAQSzOubFPbYHWE1rklgJPBQ7ZwQT6lqSxBnAqQS1kXg7CJqvHiJY10lEpAtKCj1VEXBu7Gixk9baxCrg98CDsXP5wHTaNj0toXXU0/fTErWI9HKavNabOVBL22anVcDuDsoXxMqcAgxMR4Ai0lOkOnlNNYXezICy2HFp7DEHVgB/RzChLn5dpibgtNjtMcDkuGNK7N9y9FshksP0559tDDgLmEHQD9Ey6uk6gkl1LyccDwG74p6fD0ykbcJoOUbGXl9EspaSQrZKNurp1NgRz4EG2ieLl4GnaLvz3CCSJ4vJBH0gItLrKSlkq1RGPUHwzb8kdpyVcK6ZoM9iE22TxZ+Bn9F2H+sRtDZBxR8TgX7v9iJEJN2UFKRjfQi2Ix0LXJBwrpFg69KXaZs0HieopcS/xjja9lu0HGVoTr1ID6OkIO9OP4K5EpVJzu0BNtOaKFqSxv8A++PKFRKMhErs7J5MMItb/RciaaekIN2viNYlOuI5wfLhiX0XG4BfE4yOajGM9k1RU9BwWpGIKSlI+hjBCKaRBPtLxDsKbKV9wlgO/DihbOJw2pZjPPqNFjlB+hOSniGfoBZwCvChhHMHCfapSOzw1nBakW6npCA93wCCeRczkpyLH04bnzQ0nFbkXVFSkN6tGJgbO+K1DKeNr1lsouPhtMlGR2k4reQgJQXJTvHDac9POBc/nDb++A3ww4TXaBlOm5g0NJxWspSSguSe4xlO23J0NZw2PmloOK30YkoKIvE6G077Nu07u1MdTjsZmISG00qPF1lSMLMygtX+TyZowV3q7t9JKGPAdwjGmxwErnH3F6OKSeRdM4Lf5JNJPpx2G+07u5eT+nDacoKlzUUyLMqawlHgFnd/0cwGA6vN7Cl33xBX5oME358mAe8B7o39K9J7tAyFnUjwGx2vZTht4gipZMNpJ5C8w7ur4bTaYU+6UWRJwd13EPy64u77zKwGGE1Q4W5xCfCgBzv9rDSzoWY2MvZckd4v1eG08cfTtB1OO5DkfReTgKHA19EOe9Jt0tKnYGblwEzg+YRTo4E34u7Xxh5TUpDsdzzDaV8m2FXv57QdThvv3tjRl+CrV1nstshxiDwpmNkg4BfAInffm3g6yVPa7Q9qZgsIdgZg7Nix3R6jSI/S1XDaLQTNUKuAn8TuxyeKIwQjo1qWFSknGFobf7Q8NiCaS5DeK9KkYGYFBAlhmbs/kqRILcH3mRZjgDcTC7n7UoKtYpg9e3bv2lRapDv1AypixyVAPcFfRsvfM2qCAAAIJklEQVQOe58Erifo+N5GsJ7UNoI6+sO0HSUFwT4aiYki/hiKhtfmmChHHxnBVKAad7+rg2KPAQvN7KcEHcx71J8gchyS7bB3TgdljxGsUruV1qTRcmwAfgccSnjOEDquZYwjmA2upJFVoqwpzAOuAtaa2ZrYY7cRVIpx9yXAbwmGo75CME7j2gjjEck+qe6wB5BH0GM3muCvM5ET1DwSaxktx3PA7oTnFBL8RZeTPHmMir2v9BpRjj56ji6+Q8RGHd0YVQwichwMKI0diZP3WuyhfS1ja+zfNcDOhPL5BI3CHTVRlaH1pXoYzWgWkdQV0fEQWwian14neRPVswQ9hvGd4i2d4Z01UWkWeFopKYhI9+lPMI9iSgfnmwiGlyRronqBYFhKYmd4MR13hJejzvBupqQgIulTQLBD3vgOzjcTdJYna6KqAZ4g6H2MN5iOaxnjgJNQ0jgOSgoi0nP0obUz/Kwk551gJniyjvBtBKvZJnaG96N9DSM+eYxCn4Rx9F8hIr2HEcytKAFmdVBmL8k7wrcRrGj7dkL5PFo7w8tpnzzGklOd4UoKIpJdhgDTY0cyLZ3hyZqo/gBsp/1SIomd4eUJ9wd1Y/wZpqQgIrkllc7w7SRvnloN/JJg9ni84SSvZbQkkGH0mn4NJQURkXgFBB/k5R2cbyaYGZ6siWoT8CRwIOE5g+h8OZGT6Hp71zQtka6kICJyPPoQdE6Pov0Kt9C2MzxZE9UK2u6lAUGfxVg67gwfTdqWSFdSEBHpTql0hu+j4+VEHqd9Z3i8liXSC2m/VlU3UFIQEUm3wUBV7EjmMK2d4X8BHgA2EixqOAC4FPhWNKEpKYiI9DSFtO60dwHwKsHkvUKChDGEyPoVuuraEBGRTGtZIn1l7N+3onsr1RRERHq641ki/QSppiAiIiElBRERCSkpiIhISElBRERCSgoiIhJSUhARkZC5e6ZjOC5mVkcwz+/dKAHquzGc3kDXnBt0zbnhRK55nLuXdlWo1yWFE2Fmq9x9dqbjSCddc27QNeeGdFyzmo9ERCSkpCAiIqFcSwpLMx1ABuiac4OuOTdEfs051acgIiKdy7WagoiIdCIrk4KZfcDMNpnZK2Z2a5Lz/czsodj5582sPP1Rdq8UrvmvzexFMztqZpdlIsbulsI132xmG8zsJTN7xszGZSLO7pTCNd9gZmvNbI2ZPWdmUzMRZ3fq6prjyl1mZm5mvXpEUgo/42vMrC72M15jZp/p1gDcPasOII9gS4oJQF+CfYumJpT5ArAkdvsK4KFMx52Gay4HZgAPApdlOuY0XfM5wIDY7c/nyM95SNzti4EnMh131NccKzcY+G+CHQdmZzruiH/G1wD3RBVDNtYUzgRecffX3P0I8FPgkoQylwA/it1+GDjPzCyNMXa3Lq/Z3be6+0tAcyYCjEAq1/wHdz8Yu7sSGJPmGLtbKte8N+7uQIJt5HuzVP6eIdjW/l8J9iXrzVK93shkY1IYDbwRd7829ljSMu5+FNgDFKclumikcs3Z5niv+Xrgd5FGFL2UrtnMbjSzVwk+JL+Yptii0uU1m9lMoMzdH09nYBFJ9ff6Y7Fm0YfNrKw7A8jGpJDsG3/it6VUyvQm2XY9qUj5ms3sU8Bs4P9FGlH0Urpmd/+eu08E/h74x8ijilan12xmfYBvA7ekLaJopfIz/jVQ7u4zgKdpbfXoFtmYFGqB+Mw5BnizozJmlg8UAe+kJbpopHLN2Salazaz84HbgYvdvTFNsUXleH/OPwU+EmlE0evqmgcDVcByM9sKzAEe68WdzV3+jN29Ie53+T5gVncGkI1J4QVgkpmNN7O+BB3JjyWUeQz4dOz2ZcCzHuvB6aVSueZs0+U1x5oV/p0gIezMQIzdLZVrnhR398PA5jTGF4VOr9nd97h7ibuXu3s5Qd/Rxe6+KjPhnrBUfsYj4+5eDNR0awSZ7m2PqAf/Q8DLBL34t8ce+xrBLwtAIfBz4BXgz8CETMechms+g+BbyAGgAVif6ZjTcM1PA28Da2LHY5mOOQ3X/B1gfex6/wBMy3TMUV9zQtnl9OLRRyn+jP8l9jP+S+xnXNGd768ZzSIiEsrG5iMREXmXlBRERCSkpCAiIiElBRERCSkpiIhISElBRERCSgoiIhLKz3QAIj2Bmf0TMJ9gMbJ6YDXBQokLCJYwfgW4yt0PmtkDwCGgAhgHXEswQ34u8Ly7XxN7zf3A94DzgV3AbQSL1I0FFrn7Y7G9PH5MsKIpwEJ3r472akU6ppqC5LzYOjkfA2YCHyVYPA/gEXc/w91PJVhK4Pq4pw0DzgVuIlig7NvANGC6mZ0WKzMQWO7us4B9wB3ABcClBDNUAXYCF7j76cAngH+L5CJFUqSaggicDfzK3Q8BmNmvY49XmdkdwFBgEPD7uOf82t3dzNYCb7v72thz1xNsaLQGOAI8ESu/Fmh096bYc8pjjxcA98QSyTFgcjSXKJIaJQWR5MsVAzwAfMTd/2Jm1wDvizvXskplc9ztlvstf1dN3rqOTFjO3Ztjq/NCUNN4GziVoObe2zeJkV5OzUci8BxwkZkVmtkggtVFIViWeYeZFRD0N0ShCNjh7s3AVQTbMYpkjGoKkvPc/QUze4xg1cltwCqCTuZ/Ap6PPbaWIEl0t+8DvzCzjxOseHkggvcQSZlWSRUBzGyQu+83swEEG8AvcPcXMx2XSLqppiASWGpmUwn22viREoLkKtUUREQkpI5mEREJKSmIiEhISUFEREJKCiIiElJSEBGRkJKCiIiE/j9RdkPuwk/AbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "price_ww = []; price_bb = []; gammas = [0.0001, 0.001, 0.01, 0.05, 0.1, 0.3, 0.5]\n",
    "TC.cost_b=0.01; TC.cost_s=0.01\n",
    "\n",
    "for gamma in gammas:\n",
    "    TC.gamma = gamma\n",
    "    price_ww.append(TC.price(N=400, TYPE=\"writer\")) \n",
    "    price_bb.append(TC.price(N=400, TYPE=\"buyer\"))\n",
    "\n",
    "plt.plot(gammas, price_ww, color='green', marker='o',label=\"Writer\")\n",
    "plt.plot(gammas, price_bb, color='magenta', marker='*',label=\"Buyer\")\n",
    "plt.xlabel(\"gamma\"); plt.ylabel(\"price\"); plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we have found that:\n",
    "\n",
    "- The option pricing is an increasing function of the risk aversion coefficient for the writer, and a decreasing function for the buyer.\n",
    "\n",
    "- The option pricing is an increasing function of the transaction costs for the writer, and a decreasing function for the buyer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec5.3'></a>\n",
    "### Is the drift important? \n",
    "\n",
    "As we know from the \"classical\" no-arbitrage martingale pricing theory, the option price does not depend on the stock expected value. \n",
    "\n",
    "However, this model is a utility based model i.e. a model that does not consider a risk neutral investor. \n",
    "\n",
    "We can see that in this model the option price depends on the drift. \n",
    "\n",
    "If we consider a high risk aversion coefficient, the option price is not very sensitive to the drift. If instead we choose a small value of $\\gamma$, i.e. the investor is risk neutral, the drift plays the role of the risk neutral expected return $r$ and therefore changing $\\mu$, is like changing $r$.\n",
    "\n",
    "Following Hodges-Neuberger [2], in the practical computations **it is better to set $\\mu=r$.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt4VdWd//H3lyTcIVwSriGEmxLukYhBtFa04nTqOPOro1VqFW3pOPCrWttnrJ1f+2urz2Od6uPMWHXodLQ+pZeZSjutv2krtbVTBlFB0ABBQW6GiyTIVUgIyff3xz7ZnHNykpyQnJzk5PN6nvOcffZee5+1cmB991577bXM3REREQHole4MiIhI16GgICIiIQUFEREJKSiIiEhIQUFEREIKCiIiElJQEBGRkIKCiIiEFBRERCSUne4MtFVeXp4XFRWlOxsiIt3Khg0bqt09v7V03S4oFBUVsX79+nRnQ0SkWzGzPcmkU/ORiIiEFBRERCSkoCAiIqFud08hkbq6OiorK6mpqUl3VtKqb9++FBQUkJOTk+6siEg3lRFBobKykkGDBlFUVISZpTs7aeHuHD58mMrKSiZMmJDu7IhIN5URzUc1NTUMHz68xwYEADNj+PDhPf5qSSQTrSxfSdHjRfT6Ri+KHi9iZfnKlH1XRlwpAD06IDTS30Ak86wsX8nSXy3lVN0pAPYc28PSXy0FYPHMxR3+fRlxpSAikqm+8ruvhAGh0am6U3z1pa+m5PsUFDrIvffey+OPPx5+XrRoEZ/97GfDz/fddx+PPfZYk/0uvfRSAHbv3s2PfvSj1GdURLosd+edw+/w3JvPcdcLd1HyLyW8d/y9hGn3Htubkjz0yKCQiva5Sy+9lLVr1wLQ0NBAdXU1W7ZsCbevXbuWBQsWhJ/r6+vD9XB+QaHxGCLSPR2tOcqL777IN//4TT6+8uMMf2Q4Fz5xIbf94jZWlq8kr38euX1yE+5bmFuYkjxlzD2FZKWqfW7BggXce++9AGzZsoUZM2Zw4MABjhw5Qv/+/amoqODYsWNceeWVjB49mk2bNrF161YGDhzIyZMnuf/++6moqGDOnDncdtttfOELX+D+++/n5Zdfpra2lmXLlvH5z3+el19+mW984xsxxxCRrq++oZ4tVVtYV7kufFVUVwBgGNNHTOeTxZ+krKCMsoIypuZNJatXVpM6C6B/Tn8euuqhlOQz44LCPb+5h00HNzW7fV3lOmrra2PWnao7xZ3/eSff2/C9hPvMGTWHx699POG2RmPGjCE7O5u9e/eydu1a5s+fz759+3jllVfIzc1l1qxZ9O7dm9dee43Nmzc36Tb68MMP853vfIcXXngBgBUrVpCbm8vrr79ObW0tCxYs4JprrgFo9hgi0nUc+vBQTAB4ff/rnDxzEoC8/nmUFZSxeOZiygrKuHjsxQzuMzjhcRpPVr/60lfZe2wvhbmFPHTVQym5yQwZGBRaEx8QWlvfFgsWLGDt2rWsXbuWL37xi+zbt4+1a9eSm5sb3juYN29eUpX5iy++yFtvvcXPfvYzAI4dO8b27dvp3bt30scQkc5xpv4Mbx58MwgA+4IgsPPITgCye2UzZ9Qcbp99e3gVMHHoxDb1Flw8c3HKgkC8jAsKrZ3RFz1exJ5jTQcLHJ87npdvf7ld3914X6G8vJwZM2Ywbtw4Hn30UQYPHswdd9wBwIABA5I6lrvzz//8zyxatChm/csvv5z0MUSk47k7lccrz10F7FvHhv0bwhPLMYPGML9gPneV3kVZQRlzR8+lX06/NOc6eRkXFFrz0FUPpax9bsGCBTz66KNMnDiRrKwshg0bxtGjR9myZQvf+9732Lx5c7P7Dho0iBMnToSfFy1axFNPPcXChQvJycnhnXfeYezYse3Oo4i0zam6U2zYvyHmKmD/if0A9M3uy9zRc1k+b3l4FVAwuCDNOW6fHhcUUtk+N3PmTKqrq7nlllti1p08eZK8vLwW9501axbZ2dnMnj2b22+/nbvvvpvdu3dz0UUX4e7k5+fzi1/8ot15FJHmuTs7PtgRcxXw5sE3qfegp9+koZO4sujKMADMGjmL3lm905zrjmXunu48tElpaanHT7JTUVFBcXFxmnLUtehvIZK8YzXHeH3/66yrXMcrla/wauWrHD59GICBvQcyb+w85hfMp6ygjEvGXkL+gFYnLuuyzGyDu5e2lq7HXSmISM9U31BPRXVFTI+grVVbcYIT42n507j+wuvDq4Bp+dPI6pWV5lx3PgUFEclIVR9W8eq+V8MA8Nq+1zhxJrhvN6zfMMoKyrhp+k1hl9AhfYekOcddg4KCiHR7dfV1vPn+mzFXAe8eeReALMti9qjZ3Drr1vAqYPKwyRpAshkKCiLS7ew7vi/mZvD6/eupORsMGz9q4CjmF8xn6dyllBWUUTqmlP45/dOc4+5DQUFEurTTdad548AbMV1CK49XAtA7qzdzR88NnwkoKyhj3OBxugpoBwUFEeky3J2dR3bGXAVsOriJsw1nAZgwZAKXF14eBoDZI2fTJ7tPmnOdWRQUOkhWVhYzZ87E3cnKyuKJJ54Ih7YQkcSO1x7n9X2vx1wFVJ+qBmBAzgDmjZ3Hly/9ctgldOTAkWnOcebruUHhAPAp4KfAqPYfrl+/fmzaFAzE99vf/pavfOUr/PGPf2z/gZtRX19PVlbP6y4n3VeDN7Ctelt4FfBK5StsObQl7BI6NW8qn7jgE5SNDa4Cpo+YTnavnltFpUvP/Yt/C1gDfBN4smMPffz4cYYOHQoEYxVFj366fPlySktLGTduHE888QQ///nPAVi9ejVPPfUUq1at4sUXX+TrX/86tbW1TJo0iWeeeYaBAwdSVFTEHXfcwYsvvsjy5cv51Kc+1bEZF+lAh08djukS+uq+VzleexyAIX2HUFZQxg3FN1BWUMa8sfMY2m9omnMskIlB4R6g+ZGz4U9AQ9TnpyKvXsDlzewzB2h5nD1Onz7NnDlzqKmp4cCBA/z+979vMf3ChQtZtmwZVVVV5Ofn88wzz7BkyRKqq6t58MEH+d3vfseAAQP49re/zWOPPcbXvvY1APr27cuaNWtazoxIJ6urr6P8UHlMl9DtH2wHoJf1YtbIWdwy45bwXsCU4VPoZT1yjq8uL/OCQmvmATuBaoLg0AvIAya177DRzUevvPIKn/nMZ1ocAM/MuPXWW/nhD3/IkiVLeOWVV3juuef4zW9+w9atW8NZ2s6cOcP8+fPD/W666ab2ZVSkA+w/sT8mAKzfv57TZ08DMHLASOaPm8+dJXcGo4SOmcvA3gPTnGNJVsqCgpmNA54jaLFvAFa4+z/GpckFfggURvLyHXd/pl1f3MoZPQB3ASuAvsAZ4JN0aBPS/Pnzqa6upqqqiuzsbBoazl2a1NTUhMtLlizhuuuuo2/fvvz1X/812dnZuDsf+9jH+PGPf5zw2Bo2WzrayvKVLQ4QWXO2ho0HNsbcDG6cHzinVw4Xjb4ofCagrKCM8bnj1SW0G0vllcJZ4D53f8PMBgEbzGy1u0fPH7kM2Oru15lZPvC2ma109zMpzBe8D/wNsJQgOBzo2MNv27aN+vp6hg8fzvjx49m6dSu1tbXU1NTw0ksvcdlllwHBbG1jxozhwQcfZPXq1QCUlZWxbNkyduzYweTJkzl16hSVlZVccMEFHZtJERJPT/u5X36OtXvXkt0rm3X71rHxwEbqGuqAYN6R+QXzubfsXsoKypgzag59s/umswjSwVIWFNz9AJHq1t1PmFkFMBaIDgoODLLgtGIg8AFBMEmtVVHL3+2YQzbeU4Cgr/UPfvADsrKyGDduHDfeeCOzZs1iypQplJSUxOy3ePFiqqqqmDZtGgD5+fk8++yz3HzzzdTWBpN2PPjggwoK0uHO1J/hyy9+OWZuEYDTZ0/z5Pon6Z/Tn4vHXMwX538x7BI6etDoNOVWOkun3FMwsyKgBHg1btMTwC+B/cAg4CZ3b6Abqq+vb3bbI488wiOPPJJw25o1a/jc5z4Xs27hwoW8/vrrTdLu3r27XXmUnulozVG2VW9jW/U2Kqoq2HY4eN95ZGc4T0A8wzh2/zF1Ce2BUv6Lm9lA4HngHnc/Hrd5EUFfoYUEt3pXm9mf4tOZ2VKCxh4KCwtTneVOM3fuXAYMGMCjjz6a7qxIN+fu7DuxL6j0q7dRUX3u/eDJg2G6nF45XDD8AmaNnMVN02/iqfVPhfMHRCvMLVRA6KFS+qubWQ5BQFjp7qsSJFkCPOzBTD87zGwXMBV4LTqRu68gaP2ntLS0e80K1IINGzakOwvSzdTV17Hjgx0xlX7jVcDJMyfDdLl9cinOL+bayddSnFfM1LypFOcVM2HohJjKfmr+1JRNTyvdUyp7HxnwfaDC3R9rJtle4CrgT2Y2EriQoMNom7l7j+/x0N1m0ZPmHa89nrDJ590j74bjAAEUDC6gOK+YJXOWnKv884sZOWBkUv8fUjk9rXRPKZuO08wuI3hUrJxzj4s9QND9FHd/2szGAM8CowEjuGr4YUvHTTQd565duxg0aBDDhw/vsYHB3Tl8+DAnTpxgwoQJ6c6OJMHdOXDyQMImn8aJ4QGye2UzZdiU8Gy/seK/cPiFDOozKI0lkO4k7dNxuvsagoq+pTT7gWva+10FBQVUVlZSVVXV3kN1a3379qWgoCDd2ZA4dfV17DyyM2GTT+OwDwCDeg+iOL+YqydeHdPkM3HoRHKyctJYAulJMuJOUk5Ojs6OJe1OnjmZsMlnxwc7wn7+AGMGjaE4r5hbZ90ac+Y/euDoHnulK11HRgQFkc7i7rz/4fsJm3waJ36BYArIycMmU5xfzPUXXk9xflD5T82byuA+g9NYApGWKSiIJHC24Sy7juxK2ORztOZomG5g74FMzZvKR4s+GtPkM2nYJHpn9U5jCUTOj4KC9GgfnvmQtw+/3aTJZ/sH2zlTf260lVEDR1GcV8zNM26OafIZO2ismnwkoygoSMZzd6pOVSVs8mkc2A2CIZ4nDZ1EcX4xfz7lz2OafIb0HZLGEoh0HgUFyRj1DfXsPrr7XKUfdeZ/pOZImK5/Tn+m5k3l8sLLY7p5Th42WfP9So+noCDdzum60zFNPo1B4J3D71BbXxumGzlgJFPzpnLT9JvC5p6peVMpGFygCV5EmqGgIF1W9anqhE0+e47uCef17WW9mDBkAsX5xSyatCimyWdYv2FpLoFI96OgICnV2gQuDd7AnqN7Ejb5RA/U1i+7HxfmXcj8gvkxQzpMGT5F4/mLdKCUDXORKomGuZCuKX4CF4DeWb35xJRP0Du7NxVVFbx9+G1qzp6bjS6/f36T4Rym5k2lMLdQTT4i7ZD2YS6k5zp48iCbDm5i+f9b3mQClzP1Z1i1bVXY5HPVhKsozi8Og8Dw/sPTlGsRAQUFaYcGb2DXkV1sPLiRjQc2Bu8HN8aM35+IYey8+7wGwxWRFFNQkKTU1dextWprTAB48/03wwHdsiyLafnTuGbSNZSMKqFkVAmf/vmnY4Z+aFSYmzkTJYlkGgUFaeLkmZO8efBNNh7cyKaDm9h4cCObD20On/Dtn9Of2SNn8+mZn6ZkdBAApo+Y3uSG78NXP6wJXES6GQWFHu7Qh4eCij+q+Wf74e1hl8+8/nmUjCrhnkvuYc6oOZSMLmHKsClk9cpq9diawEWk+1Hvox7C3dl9dHeT9v/oyVyKhhQFFX+k+adkdInG9hHJEOp91IPV1dexrXpbTADYdHATx2qPAUH7/9S8qSycsDAMAHNGzWFov6FpzrmIpJuCQjf34ZkPeev9t2La/8vfLw+He+iX3Y9ZI2dx84ybKRkdVP4zR8ykX06/NOdcRLoiBYUupLWnf6tPVTdp/3/n8Ds0eDAF9tC+QykZXcLyecvD5p8Lhl9Adi/9zCKSHN1T6CISPf3bJ6sPH5/8cc76WTYe3BjTvXPc4HFhz5/GADBu8Di1/4tIQrqn0E0crTnKlkNb+MKvv9Dk6d/a+lp+/vbPKc4r5iPjPxK2/ZeMKtGTvyKSEgoKnaTmbA0VVRVsPrSZ8kPl4Xuih7uiGcbWZVs7KZci0tMpKHSw+oZ63j3yLuXvBxX/5qrNlL9fzvYPtodt/72zelOcV8wV469g5oiZzBgxg8+/8Hn2ndjX5Hh6+ldEOpOCwnlyd/ad2BdU/FFn/1urtoajfhrGpGGTmDliJjdOvzEMAJOHTSYnKyfmeN/+2Lf19K+IpJ2CQhKOnD7SpNln86HNHK05GqYZM2gMM0bMYNnFy5gxYgYzRsxgWv40+uf0T+o79PSviHQFPaL3UWtdPRudrjvN1qqtTQJA9FO/uX1ymTlyJjPyZwTvI2YwPX+6bvyKSJem3kcR8V099xzbw9JfLWX/if0U5RbFBIAdH+wIx/zpk9WHafnTuHri1czID878Z46cqWEfRCSjpexKwczGAc8Bo4AGYIW7/2OCdB8FHgdygGp3v6Kl47b1SqHo8SL2HNvT7PZe1ospw6aETT6N7f6Thk3SQ18ikjG6wpXCWeA+d3/DzAYBG8xstbuH/SvNbAjwJHCtu+81sxEdnYm9x/Y2u+2NpW8wNW+qhnwQEYlI2aS37n7A3d+ILJ8AKoCxccluAVa5+95IukMdnY/munSOzx1PyegSBQQRkSidMhO6mRUBJcCrcZsuAIaa2ctmtsHMPtPM/kvNbL2Zra+qqmrTdz901UNNegCpq6eISGIpDwpmNhB4HrjH3Y/Hbc4G5gJ/DiwC/o+ZXRB/DHdf4e6l7l6an5/fpu9fPHMxK65bwfjc8RjG+NzxrLhuhbp6iogkkNI7qWaWQxAQVrr7qgRJKgluLn8IfGhm/w3MBt7pyHwsnrlYQUBEJAkpu1KwoN/m94EKd3+smWT/CVxuZtlm1h+4hODeg4iIpEEqrxQWALcC5Wa2KbLuAaAQwN2fdvcKM/sN8BZBt9V/dffNKcyTiIi0IGVBwd3XAK0+5eXu/wD8Q6ryISIiyeuU3kciItI9KCiIiEhIQUFEREIKCiIiElJQEBGRkIKCiIiEFBRERCSkoCAiIiEFBRERCSkoiIhISEFBRERCCgoiIhJSUBARkZCCgoiIhBQUREQkpKAgIiIhBQUREQkpKIiISEhBQUREQgoKIiISUlAQEZGQgoKIiIQUFEREJKSgICIioaSDgpmNN7OrI8v9zGxQ6rIlIiLpkFRQMLPPAT8D/iWyqgD4RaoyJSIi6ZHslcIyYAFwHMDdtwMjWtrBzMaZ2R/MrMLMtpjZ3S2kvdjM6s3shmQzLiIiHS87yXS17n7GzAAws2zAW9nnLHCfu78RaWraYGar3X1rdCIzywK+Dfy2bVkXEZGOluyVwh/N7AGgn5l9DPgP4Fct7eDuB9z9jcjyCaACGJsg6f8GngcOJZ1rERFJiWSDwv1AFVAOfB74L+Dvk/0SMysCSoBX49aPBf4KeDrZY4mISOok23zUD/g3d/8ehE0+/YBTre1oZgMJrgTucffjcZsfB/7O3esbm6aaOcZSYClAYWFhklkWEZG2SvZK4SWCINCoH/C71nYysxyCgLDS3VclSFIK/MTMdgM3AE+a2V/GJ3L3Fe5e6u6l+fn5SWZZRETaKtkrhb7ufrLxg7ufNLP+Le1gwan/94EKd38sURp3nxCV/lngBXdXV1cRkTRJNih8aGYXNd44NrO5wOlW9lkA3AqUm9mmyLoHgEIAd9d9BBGRLibZoHAP8B9mtj/yeTRwU0s7uPsaoPkbBU3T355sWhERSY2kgoK7v25mU4ELCSr6be5el9KciYhIp2sxKJjZQnf/vZn9r7hNU8yMZm4ei4hIN9XalcIVwO+B6xJsc0BBQUQkg7QYFNz962bWC/i1u/97J+VJRETSpNXnFNy9AVjeCXkREZE0S/bhtdVm9qXIyKfDGl8pzZmIiHS6ZLuk3kFwD+Fv49ZP7NjsiIhIOiUbFKYRBITLCILDn9AgdiIiGSfZoPADggl2/iny+ebIuhtTkSkREUmPZIPChe4+O+rzH8zszVRkSERE0ifZG80bzays8YOZXQL8T2qyJCIi6ZLslcIlwGfMbG/kcyFQYWblgLv7rJTkTkREOlWyQeHalOZCRES6hGQHxNuT6oyIiEj6JXtPQUREegAFBRERCSkoiIhISEFBRERCCgoiIhJSUBARkZCCgoiIhBQUREQkpKAgIiIhBQUREQkpKIiISEhBQUREQgoKIiISSllQMLNxZvYHM6swsy1mdneCNIvN7K3Ia62ZzU50LBER6RzJzqdwPs4C97n7G2Y2CNhgZqvdfWtUml3AFe5+xMz+DFhBMKGPiIikQcqCgrsfAA5Elk+YWQUwFtgalWZt1C7rgIJU5UdERFrXKfcUzKwIKAFebSHZncCvOyM/IiKSWCqbjwAws4HA88A97n68mTRXEgSFy5rZvhRYClBYWJiinIqISEqvFMwshyAgrHT3Vc2kmQX8K3C9ux9OlMbdV7h7qbuX5ufnpy7DIiI9XCp7HxnwfaDC3R9rJk0hsAq41d3fSVVeREQkOalsPloA3AqUm9mmyLoHgEIAd38a+BowHHgyiCGcdffSFOZJRERakMreR2sAayXNZ4HPpioPIiLSNnqiWUREQgoKIiISUlAQEZGQgoKIiIQUFEREJKSgICIiIQUFEREJKSiIiEhIQUFEREIKCiIiElJQEBGRkIKCiIiEFBRERCSkoCAiIiEFBRERCSkoiIhISEFBRERCCgoiIhJSUBARkZCCgoiIhBQUREQkpKAgIiIhBQUREQkpKIiISEhBQUREQgoKIiISUlAQEZFQyoKCmY0zsz+YWYWZbTGzuxOkMTP7JzPbYWZvmdlFqcqPiEi3dgC4AjiY2q9J5ZXCWeA+dy8GyoBlZjYtLs2fAVMir6XAUynMj4hI9/UtYA3wzdR+TcqCgrsfcPc3IssngApgbFyy64HnPLAOGGJmo1OVJxGRLu8MsBv4E/BjIAcwglPmhsi7Af1S8/XZqTlsLDMrAkqAV+M2jQXei/pcGVl3IG7/pQRXEhQWFqYqmyIiqVUH7COo6d5r5v19wOP2ywbqI+v7A38FfCc1WUx5UDCzgcDzwD3ufjx+c4Jd4v8cuPsKYAVAaWlpk+0iImlXB+yn5Qr/IE1ruEHAuMhrduS9IOq9APg7ghqwD1ADDAZGpaYYKQ0KZpZDEBBWuvuqBEkqCYreqIDgzyoi0nXUEbRfNFfZv0fiCn8g5yr8mTSt8McRVPCteR/4G4L2khXEtaV0rJQFBTMz4PtAhbs/1kyyXwLLzewnwCXAMXdPYXFFROKcpWmFn+gMvyFuv8YKvwCYQdPKvgDI7aA8Rp9Sf7eDjtmMVF4pLABuBcrNbFNk3QNAIYC7Pw38F/BxYAdwCliSwvyISE/TWOG31KRzgKYV/gDOneFPp/kz/EQN4N1cyoKCu6+hlT+ZuzuwLFV5EJEMdpbgDL6lJp1EFX5/zlX4HyNxhZ9LRlb4yeiU3kciIm1ST3IVfn3cfo0VfgFBhZ+oSWcIPbbCT4aCgoh0rsYKv6Umnf00rfD7ce4M/yoSn+Grwm83BQUR6Tj1BD1l4iv66OXmKvzGiv1KzlX+0ZX+UFThdwIFBREJmmI+BfyU5vu/NxBU+C016ewnaOuP1pdzFfuVJG7SGYYq/C5CQUFEgvF0/gQsB24hcaW/j8QVfmPlfgWJm3RU4XcrCgoiPYEDR4BdkdfOyPsKYnvnPB95QfD0bGPFfjmJm3SGowo/wygoiGSKGoKB1Bor/F1xy8fi0g8jeOjqCEHz0VmCQLAQeISgf74q/B5HQUGku6gnaMJJVOHvpOnQB32BCZHXZVHLEyPvjcMr3EVwxdCXYITOIoJgIT2SgoJIV+HABySu8HcBewjG4GnUi6AJZyJwLbEV/gSCG8bJnOl34rg60vUpKIh0plMETTzNne2fiEufR1DBXwTcQOzZ/jigdwfkqRPH1ZGuT0FBpCPVE/TUae5sP34qxX6cq+Q/QtOz/UGdkmuRkIKCSFs4cJjmb+buIbbbZi+CISAnEAz9GF3hTwRGoJu50qUoKIjEO0XzzTu7gJNx6fMJKvmLgRuJPdsfRzCdokg3oaAgPc9Zgiae5s72349L359zlfyVxJ7tTyAYV18kQygoSOZxoIrmz/bfI7aJJ4tzTTzX0bTrZj5q4pEeQ0FBuqcPab55Z1dke7QRBJV8GXAzsWf749D/BJEI/VeQ1EtmsLV4dQRn9M2d7VfFpR/IubP7q4k92y8imElLRFqloCCp9y1gDcGga09G1jlwiObP9t8jdnjlbIImnonAX9K062YeauIR6QAKCtLxThM8mTuRYNiERk9FXkbQP/9U3H6jCCr4S2nadXMs+tcq0gn030yad5agcv+AoG/+4SSXT7dwzCEE7frFxJ7tFxH08hGRtFJQ6GrOp/29NU4wQmZ85d1aJR8/qma0bIJRNodH3scTDMXQuK5x/XPACwTDMdQR3OR9MsHxRKRLUFDoahK1vzdygrPwZCv1xuUjNJ3+MNoQzlXkecCFnKvUhzezPJjk2vB/RDAKpwZbE+kWzN3TnYc2KS0t9fXr16c7G+1TBxyNeh0BPkHsCJiNehGMa99Ywde2cNwBNF+RN1fJDyXopy8iGc3MNrh7aWvpes6VQkc2yzQQjGbZWKEfbWE50br4PvTNGUzQ9j4amEfLZ+7DCMbDFxFph54TFKKbZb5LMEtVaxV6c8vHiJ3CMJ4BuQTNMkMj71Oiloc0s/wwQXNLb4JeO4tR+7uIdKrMDwr9CAJAo8Zuka3pT2yFPRqYRvMVevTyYIJmn7b6EE12IiJplbKgYGb/RtBSfsjdm0zuZ2a5wA8JHknKBr7j7s90eEZ2Al8imIy8lqD9/EKCB6AKSVy5D6FjJi9pK012IiJplsorhWeBJwg6JSayDNjq7teZWT7wtpmtdPczzaQ/P6MJztzrODcH7RXAQx36LSIiGeF8GjmS4u7/TdAxstkkwCAzM4KRaz4gduzKjtM4B+26yHv87FciIgKk957CE8Avgf0Ekw7e5O4t3b49f2qWERFJSsquFJKwCNgEjAHmAE+Y2eBECc1sqZmtN7P1VVXxw2OKiEhuG7hmAAAFiklEQVRHSWdQWAKs8sAOgrExpyZK6O4r3L3U3Uvz8/M7NZMiIj1JOoPCXuAqADMbSdAnaGca8yMi0uOlskvqj4GPAnlmVgl8ncgU5u7+NMHjZM+aWTnB415/5+7VqcqPiIi0LmVBwd1vbmX7fuCaVH2/iIi0XTqbj0REpIvpdqOkmlkVsCdudR6QSU1PmVYeyLwyZVp5IPPKlGnlgfaVaby7t9pTp9sFhUTMbH0yQ8J2F5lWHsi8MmVaeSDzypRp5YHOKZOaj0REJKSgICIioUwJCivSnYEOlmnlgcwrU6aVBzKvTJlWHuiEMmXEPQUREekYmXKlICIiHaDLBwUzu9bM3jazHWZ2f4Ltfczsp5Htr5pZUdS2r0TWv21mizoz38053/KYWZGZnTazTZHX052d90SSKM9HzOwNMztrZjfEbbvNzLZHXrd1Xq5b1s4y1Uf9Rr/svFw3L4nyfNHMtprZW2b2kpmNj9rWXX+jlsrUHX+jvzGz8kie15jZtKhtHVvPuXuXfRHMk/YuMJFgLrQ3gWlxaf4WeDqy/Cngp5HlaZH0fYAJkeNkdePyFAGb0/2bnEd5ioBZBJMt3RC1fhjBWFfDCOa92wkM7c5limw7me4ynEd5rgT6R5bvivo3151/o4Rl6sa/0eCo5b8AfhNZ7vB6rqtfKcwDdrj7Tg9mZPsJcH1cmuuBH0SWfwZcFZm453rgJ+5e6+67gB2R46VTe8rTFbVaHnff7e5vAfFzZSwCVrv7B+5+BFgNXNsZmW5Fe8rUFSVTnj+4+6nIx3VAQWS5O/9GzZWpK0qmPMejPg4gmKQMUlDPdfWgMBZ4L+pzZWRdwjTufhY4BgxPct/O1p7yAEwws41m9kczuzzVmU1Ce/7GXfH3gfbnq29k7o91ZvaXHZu189LW8twJ/Po89+0s7SkTdNPfyMyWmdm7wCPAF9qyb1ukc+a1ZCQ6Q47vLtVcmmT27WztKc8BoNDdD5vZXOAXZjY97gyis7Xnb9wVfx9of74K3X2/mU0Efm9m5e7+bgfl7XwkXR4z+zRQSjCLeZv27WTtKRN009/I3b8LfNfMbgH+Hrgt2X3boqtfKVQC46I+FxBM35kwjZllA7kE8z0ns29nO+/yRC4PDwO4+waCtsMLUp7jlrXnb9wVfx9oZ748GP0Xd98JvAyUdGTmzkNS5TGzq4GvAn/h7rVt2TcN2lOmbvsbRfkJ0HiF0/G/UbpvsrRyAyab4ObWBM7dgJkel2YZsTdm/z2yPJ3YGzA7Sf+N5vaUJ78x/wQ3pPYBw7p6eaLSPkvTG827CG5gDo0sp7U8HVCmoUCfyHIesJ24G4ZdsTwEleK7wJS49d32N2qhTN31N5oStXwdsD6y3OH1XFp/3CT/YB8H3on8wF+NrPsmQfQH6Av8B8ENlteAiVH7fjWy39vAn6W7LO0pD/BJYEvkH8AbwHXpLkuS5bmY4GzmQ+AwsCVq3zsi5dwBLEl3WdpbJuBSoDzyG5UDd6a7LEmW53fA+wRzpm8CfpkBv1HCMnXj3+gfI///NwF/ICpodHQ9pyeaRUQk1NXvKYiISCdSUBARkZCCgoiIhBQUREQkpKAgIiIhBQUREQkpKIiISEhBQaQNIvNabDOzfzWzzWa20syuNrP/icw5MM/M/q+ZfSlqn83R83yIdGUKCiJtN5ngCdNZwFTgFuAy4EvAA2nMl0i7KSiItN0udy939waCoQde8mBogHKCCXhEui0FBZG2q41aboj63EAwuNlZYv9v9e2kfIm0m4KCSMfbDVwEYGYXEYxeKdItKCiIdLzngWFmtolgfuB30pwfkaRplFQREQnpSkFEREIKCiIiElJQEBGRkIKCiIiEFBRERCSkoCAiIiEFBRERCSkoiIhI6P8DtoyKZND8VuAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "price_mu1 = []; price_mu2 = []; mus = [0.01, 0.05, 0.1, 0.2, 0.3]\n",
    "TC.gamma = 1 # high value of risk aversion\n",
    "TC.cost_b=0.01; TC.cost_s=0.01\n",
    "\n",
    "for mu in mus:\n",
    "    TC.mu = mu\n",
    "    price_mu1.append(TC.price(N=400, TYPE=\"writer\")) \n",
    "    price_mu2.append(TC.price(N=400, TYPE=\"buyer\"))\n",
    "\n",
    "plt.plot(mus, price_mu1, color='green', marker='o',label=\"Writer\")\n",
    "plt.plot(mus, price_mu2, color='magenta', marker='*',label=\"Buyer\")\n",
    "plt.xlabel(\"mu\"); plt.ylabel(\"price\"); plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other references\n",
    "\n",
    "[1] Cantarutti, N., Guerra, J., Guerra, M., and Grossinho, M. (2019). Option pricing in exponential Lévy models with transaction costs. [*ArXiv*](https://arxiv.org/abs/1611.00389). \n",
    "\n",
    "[2] Hodges, S. D. and Neuberger, A. (1989). Optimal replication of contingent claims under transaction costs. The Review of Future Markets, 8(2):222–239."
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}