2.3 American Options.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# American options\n",
    "\n",
    "\n",
    "## Contents\n",
    "   - [Binomial tree](#sec1)\n",
    "   - [Longstaff - Schwartz](#sec2)\n",
    "   \n",
    "   \n",
    "The American option problem is an optimal stopping problem formulated as follows:\n",
    "\n",
    "$$ V(t,s) =  \\sup_{\\tau \\in [t,T]} \\mathbb{E}^{\\mathbb{Q}}\\biggl[ e^{-r(T-\\tau)} \\Phi(S_{\\tau}) \\bigg| S_t=s \\biggr]$$\n",
    "\n",
    "where $\\Phi(\\cdot)$ is the payoff, also called *intrinsic value*.    \n",
    "This formula resembles the formula of a European option, except that in this case the option can be exercised at any time. But... at which time? The problem is to find the best exercise time $\\tau$ that maximizes the expectation.\n",
    "\n",
    "Notice that the time $\\tau$ is a random variable! It can be different for different paths of the stock process.   \n",
    "The solution of the problem provides a strategy, that at each time assesses if it is optimal to exercise the option or not, depending on the current value of the underlying stock.\n",
    "\n",
    "This notebook is dedicated only to the pricing problem, i.e. just to find $V(t,s)$, and not to find $\\tau$.\n",
    "\n",
    "It can be proved that the function $V(t,s)$ is the (weak) solution of the following PDE:\n",
    "\n",
    "\\begin{equation}\n",
    "\\max \\biggl\\{  \\frac{\\partial  V(t,s)}{\\partial t}  \n",
    "          + r\\,s \\frac{\\partial V(t,s)}{\\partial s}\n",
    "          + \\frac{1}{2} \\sigma^2 s^2 \\frac{\\partial^2  V(t,s)}{\\partial s^2} - r  V(t,s), \\;  \n",
    "          \\Phi(s) - V(t,s) \\biggr\\}  = 0.\n",
    "\\end{equation}\n",
    "\n",
    "$$ V(T,s) = \\Phi(s) $$\n",
    "\n",
    "The numerical algorithm for solving this PDE is almost identical to the algorithm proposed in the notebook **2.1** for the Black-Scholes PDE. \n",
    "\n",
    "Let us have a look at the differences between the implementation of the European and the American algorithm, in the class `BS_pricer`:\n",
    "\n",
    "```python\n",
    "if self.exercise==\"European\":        \n",
    "    for i in range(Ntime-2,-1,-1):\n",
    "        offset[0] = a * V[0,i]  \n",
    "        offset[-1] = c * V[-1,i]\n",
    "        V[1:-1,i] = spsolve( D, (V[1:-1,i+1] - offset) )\n",
    "elif self.exercise==\"American\":\n",
    "    for i in range(Ntime-2,-1,-1):\n",
    "        offset[0] = a * V[0,i]\n",
    "        offset[-1] = c * V[-1,i]\n",
    "        V[1:-1,i] = np.maximum( spsolve( D, (V[1:-1,i+1] - offset) ), Payoff[1:-1])\n",
    "```\n",
    "\n",
    "The European and the American algorithms just differ in one line!!    \n",
    "\n",
    "I don't consider dividends in the notebook. In absence of dividends the American call has the same value of the European call. For this reason, I'm going to consider only the pricing problem of a put option\n",
    "\n",
    "$$ \\Phi(S_{\\tau}) = ( K - S_{\\tau} )^+ $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functions.BS_pricer import BS_pricer\n",
    "from functions.Parameters import Option_param\n",
    "from functions.Processes import Diffusion_process\n",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as ss\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from IPython.display import display\n",
    "import sympy; sympy.init_printing()\n",
    "\n",
    "def display_matrix(m):\n",
    "    display(sympy.Matrix(m))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creates the object with the parameters of the option\n",
    "opt_param = Option_param(S0=100, K=100, T=1, exercise=\"American\", payoff=\"put\" )\n",
    "# Creates the object with the parameters of the process\n",
    "diff_param = Diffusion_process(r=0.1, sig=0.2)\n",
    "# Creates the pricer object\n",
    "BS = BS_pricer(opt_param, diff_param)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The class `BS_pricer` implements two algorithms to price American options:\n",
    "- the Longstaff-Schwartz Method (see section [below](#sec2))\n",
    "- the PDE method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAK0AAAASCAYAAAApM17jAAAABHNCSVQICAgIfAhkiAAABrFJREFUaIHt2mmsXVUVB/Bf66vSagEBLVEJIJhQ2gRErSjBDgwKWFIw+IHIEAYnIsWpTBpunBAkKIgDDTgQookDlVAGQVIsFAiioFRQwfYpBJpKtVBLaW37/LD28Z537jnn3nPf6xdz/8nN6Vtr7el/1ll77bXLAAP8n+EUjKTfWX20Pw534hlswir8FO8qsb0Md+PpZPtPPIJLsPs4zXf3JF+Cp9I4L+A+nImJJW2Gc30Wf2sq5jEBZ+BBbMBLaS3n4hUl9qfXjJH9tlWM1YRjeBO+h2exOa3vG3hthX2Gw/FzPJfaPZfGPbZg1w/HjdYyoaaDvfCYIPk1OBvXdVlYHpdhEdbhF3ge++N4DOFU3Jiz34Lf4XGsxatxKN4uCD5UOPRY5vtRfEcQvgx/xzSciF3ESzlJOEmGYewqXmwR/8YVJfIbxAe0FrdgI47EgRVjHIwFFes6HPNwK95f0DXleD/cj9fjZvwJszAXf8Zhqa8iPocvpv6XCv72wFsFj4tytv1w3M9aOjABv8Jf8TXNI+2eIjKsEQTlMTf1t6og36miry8n+2+Pw3znYb7Or31PQe4IPlDQDadfr1igvb49cvJJIvqMiMjaKx5IbY4vyPvh+JdJ/omC/Mok/27J+JmD3YWpJfpJhb/74biftXRgIbbjPWhp7rTvTG1urtC/KLbNXnCQNmlVGOt84aLU7psF+bBmTntD6uecEt3MpPttj31l9s/oTCuacvzmZL9ap0NNFbvGRrHDZZgonGUjXtfjnOtQxXGjtZTlF9PxVVyF5X1O7kmx3c8yOtoQjjVVRMZeMD89/1ChH4/5wn/Sc2uJ7lX4kCB9ofj6y3JTImpQHhky2SEi5eiGj6Tn9Tpz2qYcz0vPO8UHnscGrMAUkYZleDf2xW34l8g5zxccVOXMdajieEz+MoSHRX4zOcla+otc5wly1mIxLsVP8LIgrrgNZPhMGvPruDeN/XvlX/p4zXdI5MMjeG9BN6z8YLQKs0v6+lHSf7xEl0XOEaOdowyThaNsE/l6GZpwnKVNn67o65qk/1hO9skku0YEjSIHv9Z7BK7juOlaRuELgqT8V9TSf/VggagC5Bf6JE6uabOmYH+7SOR35HyvSG1uLdFdIqLUNBGJZorcb7uoChxUsD859fUUdsvJh8QhJFvXMV3mdFqyW9rFrleOF6vnJTs7XJiTXZpkW1OfR4hD7gzckXT3dJlfhjqOm67lf5iVJnd5Qd7Sn9MuSv1dKfKpKWJbzA4DxXGKmIYTRBR9NrXdEfM9N9k/YbSTdUP2EpYU5BPFdpqVxBaLysNKUcb5S9Id3aX/Fclufo1NE467Oe1Xkv6CnOxy7XJb8eOcLKo5I7qnCr1w3NhfhoRzPC7ytzxamjvtnNTmphLdFHGw2JYm1w17i7rgyh0w33OS7R+1c9FesX9qW1YiGhLb8KPCUV8Ukelt2tWAg2v6PjDZPK06d56jGcf9pAcXake7MlyX9Asr9PTG8Rx9+Muuuhe3s19ZvbKILAoVSysZblJe+qjCI8k+S9LHY77nJf1javKlGuyc2r/coM1k4cQv6SwV5XFV6rtVY9OU47PS39dW2GcR7Yic7MQk+01Fm+xDuKBC3yvHjdYylISbxQm1DIeIIvJ9Iro9UDN4hiz6VSXpmXxLD33BG9IzO0GPdb7ni4rDozhKFLKbItsSu9YPczhF1KN/qH2SLmKnZLdd9RppzvGy9DxapDD5CsJUcbGwSdziZVgutuy34JU639fM9BwuGb8Jx+PtL7Xb7X44QGfU+KB2TvfGgu4YQdgm7evZA5RvHRO1DwgrxmG+8Pmkf1j3HHZGhc3eYsscEWWwInYukb1DHDI2qE+LsqvoW7rMrSnH9He5cGPSfakgPyqNsV5n+a4JxzRcy5Cx4W7xAvc1+mv7mairHSmS7yVpQtPFVeQEsaVk+eD7xFazXNxqrRMHsdniBa8R17JjxWnaFYd7xQGhiGH8IP37pDTPZaIov0F8qMeJiHib8mvcuwTJK1ObGeKOfrPYcuui84fTc3GXtTTlmCjD3Y+rRRrwhCjszxUHxItLxvlUsrlY1EwfEu/8BMHj2cJxMzTluN+11KKlOnINJ90+JbpJIqd5UBxEtooa3FKdJ+eZ+JbYSp5Pti+IXKql2am+br6Zru53T85+Nn4s7ujXiy39H8IpT1X9fzc+K2691gtHXS2i2D5d5j5d9wNYHk04zrAXvi/+b8AW/E3k0HUc7yai8erUZp24vSqrNbc043gsaxlggAEGGGCAAQaA/wIaTbIUKTW09QAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle 4.834259780628$"
      ],
      "text/plain": [
       "4.834259780628"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BS.LSM(N=10000, paths=10000, order=2)  # Longstaff Schwartz Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANMAAAASCAYAAADBs+vIAAAABHNCSVQICAgIfAhkiAAABlhJREFUaIHtmmuIVVUUx3/aWI6ViFkORampUDlgGVkGOqeyyMrSzD5ImkQZFGX2ECvB2wvTpJIehBRkSFGZGkmaJpUliT0sGXpZzgWlhmmmfL91+rDW4Z7Zs/d57HvvDMT9w+HM7LXW3nv919mvtS9UUEEFHYLJQKs+d3rYXw+sAXYAB4BtwPvACIvuLcBLwJfAbm1zSUL9+Uj/zKcxwXYk8AHwF3BI32uA6wy9ecA6YLv68A+wGZgDnOaouwtwB7AR2APsV5v7gRMs+lNj/AifYxY7H86iSBPffEyfXBz72GTlrCN8mRqj3y4uVTGdOxsJ1F7glBTOmJgHzARagBVAMzAIuAmYAEyhbeBnA0O1vR3AeSnb2QW8aCnfG2MzG3hK+7QSGUh9gIuAAPg4ojsD+B5YCzQBJwOXATlgmv693ah/MRLcJuBdYB8wGlgIjAImIoEI8QPwhKOvI4ErgVUOP3w4g2zx9eE4q01WzqIoly++cWmDLsCnwB/Ac2RfmWqQEdsInGHIrtD6tlnKB2vbAelXpnyGfkEhKGuBUy3ybsb/3R31PKP1vGqUj6PgXx+j3uUqm5qhv1+rzY0WmQ9nkC2+ebJznNWmGM7K7YsLcXFpg+nAcWRGyJF9MF2qNh865LuRpdyFgPIMpq5IwPYBp2ews2EohUEZxVtafq/FplZl36VsI9TfQfJWJyD9YMoS3zzlH0zFcFZuX2ywxsW2zTsfeBZZXtcjS1lWbAUOA8ORmaY5IhuFrAgrPOq14STgNuAcZJBsQfptO2NcDgwAlgL/Ime6WuAgsAmZbdJirL63GOU1+jZX3mjZMKAXsDOhjbv1/QZ2f3zgE98sHPvY+HLWUb6YSBWXKuBb4FegWsty+CUgHkBmjCZgETAXeA/5cNfQfvsXRUBxCYhtQJ1Ff4bKX0ZINO2+wL1iPYxw8QJy4G8FfrTov62yeyx1hDNaK3LWikM1MuCPIWeCJAQkc+YT3zzZOPax8eGso3wxkTouT6pSNNuW1ME4jEOyX9GObwUmJdgFpBtMc5DZqC/QAyH+NWQQ70e2YlHM1XqPaj+uQg6sQ4DVKvvc0Vaj4ccqbdfEJJX/DvSOlFch2cPQfkyCb7er3soEvRAByZz5xDcrxz42Ppx1lC8mUsVlOPKRzTfKkzrowkyt73ngXKTzw4BPtD6znSgCsqd5o1ig9suN8vkU0pkmcdVIVq4Ve+o+RF9gPDIj/on4FEVXJBsYplsXIdmjeiS1/pvKrknwYYPqjU3QCxEQz1mp4+vi2McmK2ed6UtiXKqQj+MnZE9ZbAcDtVlmkfVADm7HkEEWZ+87mAapfYtR/iiF1dGG11U+PUUb/ZD7qXqLrAp4CEmtHkASLquBiylkgS6MqfsC1dlOujsWiOes1PEFN8e+Nmk560xfUsWlF/EXU9HHlqM3EY70+xzyZSqf4JAHFDeYeqr9QaP8Zi3/xmEXplZnpWxns+r3SVJUVCMfyn7ap+CjWKj15lLWC/GclTq+4Oa41DYmZ53pS2xcwmzeISQzYcMw5DLzK2RGSJPxCmcM12E+LD+coi4fhNs0Mzu0HtkeDAZOtLRfq+98ynbO1HfaTNBk5N5qMXDEodNd9Y7jjklWlDq+4Oa41DYmZ53lS0niksO9dA5Ebt3NWfZWCvvfswzZGO3QAdw/xwlIXpmG0PawGqIfso1rBR6zyJeo7Gmj/Grt105k9gPxrYb26Erh0naDRd7TUnYJkozZg3t7C4WfxXwUo2NDgN9qnsMdXx+OfeNSDGchcpTWlygS4xL3c6I0WKedGUDb2Xwpcis9GvgZOdg1IvcCNyC31rNou0cdpw8UPuARwJv6dzOSng4xUev4DGhACB+I3B11Rw60Cyx9fhC5VH4cufPapD6MR1aYuyjcZVyLbP3WIzfsLUgCog4JbqPqm1iLTBb12q8hyG/+DiFbzbgZcJq+F8XohMjKWVb4cOwbl2I4K5cvUWSJixM53KM9r7L+Flk35K5pI3KYPIrcOa3EnskK23E9eUO/DngH+AX5+I8AfyNBmYIMWBd6I1nGBmSr14L8WsO8+6kFXkEOxc3qwy7kzJXDPtMBPILc2O9EPoYGJAXbP6ZPIJNNlsRDjmycxdVhi68Px75x8eWsnL6EyBqXCiqooIIKKqigggr+P/gPgsQIrNiAnyYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle 4.815639714559457$"
      ],
      "text/plain": [
       "4.815639714559457"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BS.PDE_price((7000,5000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hVVb7/8fc3CSEhgHREQQGV3oTQRAQJKCIC6gB2rIxzdRR1dFDn5zjq2C7Xq456EUXAgmJFBOmIFAEpBqRKE0FQiiItlIT1+2NtMGACIeRknySf1/Ps5+x6zicHzTdrr73XNuccIiIiR4sJO4CIiEQnFQgREcmSCoSIiGRJBUJERLKkAiEiIllSgRARkSypQEhozGyomT1xku/R3sw25FWmbD7jRjObkctjHzWzt/M6Uy5yLDGz9mHnkIJFBUIixsy+N7M0M9tlZr+a2RgzqxZinngz+x8z2xBkWmtm/xtWnvzknKvvnJsadg4pWFQgJNIuc86VBKoAPwP/CTHLg0Ay0AIoBVwIfBNinogzs7iwM0jBpQIh+cI5txf4EKiX1XYzK2tmo81sS9DaGG1mVTNtL2dmQ8xsY7B9ZDbvc5eZLc18bCbNgU+ccxud971z7s1Mx1Yzs4+DDNvM7KWj3ntA8NlrzeySTOtPM7NRZvaLma0ys9uy+x7MrJWZfWVm281sYebTPsGprDVmtjP4jGuzeY9HzexDMxsR7LvAzBpn2v69mf3dzBYBu80sLljXMdgea2YPmdnq4Pj5h1p2ZlbHzCYGP8sKM+uV3c8ihZ8KhOQLMysB9AZmZ7NLDDAEOBM4A0gDMv+CfgsoAdQHKgF/ODVkZv8PuBFo55zLql9iNnCvmf2XmTU0M8t0bCwwGlgHVAdOB97LdGxLYAVQAXgWGJzp+HeBDcBpwJ+AJ80sJYt8pwNjgCeAcsDfgI/MrKKZJQEvApc450oB5wGpWX1Rge7AB8H7DAdGmlmxTNuvBi4Fyjjn0o869t5gexegNHAzsCfIMDF4v0rBPq+YWf1j5JDCzDmnSVNEJuB7YBewHUgHNgINM20fCjyRzbFNgF+D+SrAQaBsFvu1B34EngNmAKccI08scAcwE9gX5OkTbGsNbAHisjjuRmBVpuUSgANOBaoBGUCpTNufAoYG848CbwfzfwfeOuq9xwN9gKTge7oSSDzO9/ooMDvTcgywCWib6Xu/OYt/i47B/Aqgexbv2xuYftS6V4F/hv3fkqZwJrUgJNJ6OOfKAMWBO4EvzezUo3cysxJm9qqZrTOzHcA0oEzwl3014Bfn3K/ZfEYZoC/wlHPut+yCOOcynHMvO+faBMf8G3jDzOoGn7HO/fGv7UN+yvQ+e4LZkvhWwy/OuZ2Z9l2Hb4Ec7UygZ3B6abuZbQfOB6o453bjf0HfDmwKOvTrZPezAOsz5TnI7y2YP2zPQjVgdTb5Wh6V71p8IZQiSAVC8kXwy/lj/F/b52exy31AbaClc640cEGw3vC/7MqZWZls3v5XoCswxMza5DBPmnPu5eDYesFnnJGLTt2NQbZSmdadgW/VHG09vgVRJtOU5Jx7Osg03jnXCd9iWg68dozPPXw1mJnFAFWDLId/xGMcux44K5v1Xx6Vr6Rz7i/HeC8pxFQgJF+Y1x0oCyzLYpdS+H6H7WZWDvjnoQ3OuU3AWPz58LJmVszMLsh8sPOXcF4LfGJmLbPJ0M/8fROJQcdtn+BzvwG+xp+medrMkswsISfFxjm3HvgKeCo4phFwC/BOFru/DVxmZhcHHcUJQZ6qZlbZzLoF/QD78KfmMo7x0c3M7IqgoPULjsmuf+dorwOPm9k5wb9LIzMrj++DqWVm1wffcTEzax60sKQIUoGQSPvMzHYBO/CndPo455Zksd/zQCKwFf+LbtxR268HDuD/st6M/6V4BOfcROAmYJSZNcviM9KA/8GfLtqK74+40jm3xjmXAVwGnA38gD9l0zuHP+PV+I7tjcAn+HP2E7PItx7fufwQvr9jPXA//v/DGHwraiPwC9AO+K9jfOanQb5f8d/NFc65AznM+xzwPjAB/+8yGN/vsRO4CLgqyPET8Az+9KAUQeacHhgkUpCY2aPA2c6568LOIoWbWhAiIpKliN5laWbfAzvx51LTnXPJwfnlEfgm+fdAr2NcnSIiIiGJ6CmmoEAkO+e2Zlr3LP6ywKfNrD/+2va/RyyEiIjkShinmLoDw4L5YUCPEDKIiMhxRLoFsRZ/lYUDXnXODTKz7cGNU4f2+dU5VzaLY/vib34iKSmpWZ06x7pnSEREjjZ//vytzrmKuT0+0iM9tnHObTSzSsBEM1ue0wOdc4OAQQDJyclu3rx5kcooIlIomdm6kzk+oqeYnHMbg9fN+OvDWwA/m1kVgOB1cyQziIhI7kSsQAR3o5Y6NI+/AWcxMAo/OBnB66eRyiAiIrkXyVNMlfHDHhz6nOHOuXFmNhd438xuwd+x2jOCGUREJJciViCcc2uAxlms3wb8Yaz8E3XgwAE2bNjA3r17T/atipyEhASqVq1KsWLFjr+ziBRZBfZxhBs2bKBUqVJUr16dTM99keNwzrFt2zY2bNhAjRo1wo4jIlGswA61sXfvXsqXL6/icILMjPLly6vlJSLHVWALBKDikEv63kQkJwp0gRARkchRgTgJsbGxNGnShAYNGtCzZ0/27Nlz/INOwIsvvkjdunW59tpr2bdvHx07dqRJkyaMGDEiTz9HRCQrKhAnITExkdTUVBYvXkx8fDwDBw7M0/d/5ZVX+Pzzz3nnnXf45ptvOHDgAKmpqfTundPn2IiI5J4KRB5p27Ytq1atAqBHjx40a9aM+vXrM2jQIAAGDx7MPffcc3j/1157jXvvvReA5557jgYNGtCgQQOef/55AG6//XbWrFlDt27deOaZZ7juuutITU2lSZMmrF6d1fPmRUTyVoG9zDWzfv0gNTVv37NJEwh+Vx9Xeno6Y8eOpXPnzgC88cYblCtXjrS0NJo3b86VV17JVVddRaNGjXj22WcpVqwYQ4YM4dVXX2X+/PkMGTKEOXPm4JyjZcuWtGvXjoEDBzJu3Di++OILKlSoQMuWLRkwYACjR4/O2x9URCQbakGchLS0NJo0aUJycjJnnHEGt9xyC+D7Dho3bkyrVq1Yv349K1euJCkpiQ4dOjB69GiWL1/OgQMHaNiwITNmzODyyy8nKSmJkiVLcsUVVzB9+vSQfzIRkULSgsjpX/p57VAfRGZTp05l0qRJzJo1ixIlStC+ffvD9xzceuutPPnkk9SpU4ebbroJ8DeuiYhEI7Ug8thvv/1G2bJlKVGiBMuXL2f27NmHt7Vs2ZL169czfPhwrr76agAuuOACRo4cyZ49e9i9ezeffPIJbdu2DSu+iMhhhaIFEU06d+7MwIEDadSoEbVr16ZVq1ZHbO/VqxepqamULeufkdS0aVNuvPFGWrRoAfhWxrnnnpvvuUVEjhbRJ8rllaweGLRs2TLq1q0bUqLc69q1K/fccw8pKSc9XuFJKajfn4jknJnNd84l5/Z4nWLKJ9u3b6dWrVokJiaGXhxERHJCp5jySZkyZfjuu+/CjiEikmNqQYiISJZUIEREJEsqECIikiUVCBERyZIKxEk4NNx348aNadq0KV999RUABw8e5K677qJBgwY0bNiQ5s2bs3bt2pDTioicGF3FdBIyD7Uxfvx4HnzwQb788ktGjBjBxo0bWbRoETExMWzYsIGkpKQ8/eyMjAxiY2Pz9D1FRDJTCyKP7Nix4/Dd0Zs2baJKlSrExPivt2rVqoe3ZTZ37lzOO+88GjduTIsWLdi5cydDhw7lzjvvPLxP165dmTp1KgAlS5bkkUceoWXLljz55JP06tXr8H5Tp07lsssuA2DChAm0bt2apk2b0rNnT3bt2hWpH1tECrHC0YIIabzvQ6O57t27l02bNjFlyhTAD6dx/vnnM336dFJSUrjuuuv+MHzG/v376d27NyNGjKB58+bs2LGDxMTEY37e7t27adCgAY899hjp6enUrFmT3bt3k5SUxIgRI+jduzdbt27liSeeYNKkSSQlJfHMM8/w3HPP8cgjj5zc9yEiRY5aECfh0Cmm5cuXM27cOG644Qacc1StWpUVK1bw1FNPERMTQ0pKCpMnTz7i2BUrVlClShWaN28OQOnSpYmLO3a9jo2N5corrwQgLi6Ozp0789lnn5Gens6YMWPo3r07s2fPZunSpbRp04YmTZowbNgw1q1bF5kvQEQKtcLRgghrvO9MWrduzdatW9myZQuVKlWiePHiXHLJJVxyySVUrlyZkSNHHjHEhnMOM/vD+8TFxXHw4MHDy4eGCgdISEg4ot+hd+/evPzyy5QrV47mzZtTqlQpnHN06tSJd999N0I/qYgUFWpB5JHly5eTkZFB+fLlWbBgARs3bgT8FU2LFi3izDPPPGL/OnXqsHHjRubOnQvAzp07SU9Pp3r16qSmpnLw4EHWr1/P119/ne1ntm/fngULFvDaa68dfk51q1atmDlz5uHHn+7Zs0dDfIhIrhSOFkRIDvVBgG8RDBs2jNjYWDZv3sxtt93Gvn37AGjRosURHc8A8fHxjBgxgr/+9a+kpaWRmJjIpEmTaNOmDTVq1KBhw4Y0aNCApk2bZvv5sbGxdO3alaFDhzJs2DAAKlasyNChQ7n66qsPf/4TTzxBrVq1IvEViEghpuG+iyh9fyKFn4b7FhGRiFCBEBGRLBXoAlEQTo9FI31vIpITBbZAJCQksG3bNv2yO0HOObZt20ZCQkLYUUQkyhXYq5iqVq3Khg0b2LJlS9hRCpyEhASqVq0adgwRiXIFtkAUK1aMGjVqhB1DRKTQKrCnmEREJLIiXiDMLNbMvjGz0cFyOTObaGYrg9c/DnMqIiKhy48WxN3AskzL/YHJzrlzgMnB8jGpI1pEJP9FtECYWVXgUuD1TKu7A8OC+WFAj+O9z+7F35OxPyPvA4qISLYi3YJ4HngAOJhpXWXn3CaA4LVSVgeaWV8zm2dm80ru/4VZDW7jYPrBrHYVEZEIiFiBMLOuwGbn3PzcHO+cG+ScS3bOJe8sWYXzVw5hxrl/xR3U6SYRkfwQyRZEG6CbmX0PvAd0MLO3gZ/NrApA8Lr5eG9UqvZpTG1+PxcsfoVpLe9XkRARyQcRKxDOuQedc1Wdc9WBq4ApzrnrgFFAn2C3PsCnOXm/drOf4cuGd9Ju3v/wZft/RiSziIj8Loz7IJ4GOpnZSqBTsHxcFmO0XfAC02rdSvvpjzO181MRDSkiUtTly53UzrmpwNRgfhuQcqz9sxMTF0Obbwcyo3Ya7cc/xJeXJ9Luk355F1RERA4rcHdSx8bH0mrZUGadfiXtRt7D9GsHhh1JRKRQKnAFAiAuIY5my4fzdaWutB3+F2b0HXb8g0RE5IQUyAIBEF8ynkYrPmB+uU60fu1mZvUbEXYkEZFCpcAWCICEMgnUWT6Sb0ufT/MXruXrh0aGHUlEpNAo0AUCIKliCWouHc3ypGQaP9Wb+f8eF3YkEZFCocAXCIDSp5ei6uJxrEmsT71/XE7qc1PCjiQiUuAVigIBUKZ6GSoumMD64mdz9n3dWPzqzLAjiYgUaIWmQABUqFOBU+ZMZHOx06l2exeWvTUv7EgiIgVWoSoQAJUbn0rCjMn8FleeU/tcxMoPF4YdSUSkQCp0BQLgtBZVYfIU0mKSKNurE2vGLDv+QSIicoRCWSAAzrigOnvHTCHDYinRLYV1k1eFHUlEpEAptAUCoObF57Djo0nEuQPEXpzCj1+tCzuSiEiBUagLBMA5Peqz5Z2JlMzYQUa7Dvy84MewI4mIFAiFvkAA1L26CetfH0+Z9C3sat2RrUuP+4wiEZEir0gUCICGt7RgzYtjOHX/D/zSrCPbV28LO5KISFQrMgUCoMlf27Ls6VGcsfc7NjW6mB3rfws7kohI1CpSBQIg+e8pLHzkY87as4h19buw++ddYUcSEYlKRa5AALT8Vxfm3fcedXfOYWWdy0jbtifsSCIiUadIFgiA8wZcwazb36TR9i9ZWudy9u/cF3YkEZGoUmQLBEDb/7uG6Te8TrOtE0it3Yv0tANhRxIRiRpFukAAtBt2M1P/9BItNo1iXp3ryNiXHnYkEZGoUOQLBED7D+5gSpcBtPrhfebUv5mD6QfDjiQiEjoViECHMfcxuf3jnLf6Lb5q/BfcQRd2JBGRUKlAZNJhyj+Y3Oohzl86iBnN71GREJEiTQUiEzPoMPMJvmjcj7YLXmDGBQ+BU5EQkaJJBeIoFmO0m/8cU+vcTtuZTzPtoifCjiQiEgoViCzExBptF73MtBp9uGDSI0zvPiDsSCIi+U4FIhuxxWI4b9lgZlbtTdtR9zPz6pfCjiQikq9UII4hrngszZe/xazKPWjz3l/56pbBYUcSEck3KhDHEZ9UjHNXvMfX5TvT6o3bmP3Xd8KOJCKSL1QgciDhlOLUX/4xqae0J/mlPszt/1HYkUREIk4FIoeSKiRy9tJRLCnZiibPXMX8x8aEHUlEJKJUIE5A6dNKcsbiMXyX2IT6/7yShf8zKexIIiIRowJxgsqeeQqVU8ezrnhtzvlbNxa/Mi3sSCIiERGxAmFmCWb2tZktNLMlZvavYH05M5toZiuD17KRyhApFWqVo8zciWwsVp0z77iUZcPmhB1JRCTPRbIFsQ/o4JxrDDQBOptZK6A/MNk5dw4wOVgucCo3rETizElsi6tMlZs6s/L9b8KOJCKSpyJWIJx36IHPxYLJAd2BYcH6YUCPSGWItNObn0bMlMnsiilNuas6sWbU4rAjiYjkmYj2QZhZrJmlApuBic65OUBl59wmgOC1UjbH9jWzeWY2b8uWLZGMeVLOaHsm+8dOYb8Vp+TlHflh0ndhRxIRyRMRLRDOuQznXBOgKtDCzBqcwLGDnHPJzrnkihUrRi5kHqjZ6Sx2fjIZnCOucwo/zlgbdiQRkZOWL1cxOee2A1OBzsDPZlYFIHjdnB8ZIq1Wtzpse3ciCQf3cPDCDvw8b33YkURETkokr2KqaGZlgvlEoCOwHBgF9Al26wN8GqkM+a1u70b8+MYESqf/QlqbFLYu/insSCIiuRbJFkQV4AszWwTMxfdBjAaeBjqZ2UqgU7BcaDS8sRlrXhpLhf0b2d68I9tXbQ07kohIrsRF6o2dc4uAc7NYvw1IidTnRoNz7ziPuXtG0+CBS/ihcSdilk6h9JkF7nYPESnidCd1hDS/vz0LHx1J9T1L2dCgM7s37Qg7kojICVGBiKBW/7yYuQ98wDm7FrC6XlfStu4OO5KISI6pQETY+c90Y9Yd71B/+0yW1+nB/h17w44kIpIjKhD54IKXejHtxiE03jaZRbX/RPqe/WFHEhE5LhWIfHLhkBv4ovdAkn8aw4I6V5OxLz3sSCIix6QCkY9S3uvLpK7P02L9x8yt14eDBzLCjiQikq0cFwgzO9PMOgbziWZWKnKxCq+On93NxA5P0WrNcGY37ovLOBh2JBGRLOWoQJjZbcCHwKvBqqrAyEiFKuw6TurPxNaPcN6yN/gq+S7cQRd2JBGRP8hpC+IOoA2wA8A5t5JsRmGV4zODjjMeZdK599Mm9WW+Ov8BcCoSIhJdclog9jnnDl96Y2Zx+Gc7SC5ZjNFh7jNMrncnbWYNYGbHf4YdSUTkCDktEF+a2UNAopl1Aj4APotcrKIhJtZon/oCX5x1K22mPM7Mrk+FHUlE5LCcFoj+wBbgW+DPwOfAPyIVqiiJLRZD2yUDmVbtWtqMeYivej0fdiQRESDng/UlAm84514D/6S4YN2eSAUrSuKKx9Jq+VBmnrWXNh/cw6ybEmk95M9hxxKRIi6nLYjJ+IJwSCIwKe/jFF3xJeJotmI4syp0pfXQ25lzx7DjHyQiEkE5LRAJzrldhxaC+RKRiVR0JZSOp9HyD5hbphPJr9zMvPtHhB1JRIqwnBaI3WbW9NCCmTUD0iITqWhLKp9A7WUjWVjqfJoMuJYFjxaaB+6JSAGT0wLRD/jAzKab2XRgBHBn5GIVbaVPLUGNxaNZWiKZ+v/qxaJnx4UdSUSKoBx1Ujvn5ppZHaA2YMBy59yBiCYr4sqeUYqMheNY3bAD5/z9cpaU+Jz6d14YdiwRKUKO2YIwsw7B6xXAZUAt4BzgsmCdRFCFs8tQbu4E1sefxZl/vYzlg2eGHUlEipDjtSDaAVPwxeFoDvg4zxPJEU5tUIGMrybxc6t2nHZbF1YlTebsq5LDjiUiRcAxC4Rz7p9mFgOMdc69n0+Z5CinNzuVdVMn82u7Cyh/zUWsTZxKje6Nwo4lIoXccTupnXMHUYd06M5sU5WM8ZPZY0mUuqIjP4xfFnYkESnkcnoV00Qz+5uZVTOzcoemiCaTP6iZUoNdIyeT7mIpfmkKG6etCjuSiBRiOS0QNwP/BXwJzMs0ST6rfVktfhkxidiDB3ApKWyeuy7sSCJSSOW0QNQDXgYWAqnAf4D6kQolx1avZ31+fGMCJdJ3sLdNCtsW/Rh2JBEphHJaIIYBdYEX8cWhbrBOQtL4xnNZ+8o4yhzYzG8tOrL9u81hRxKRQiano7nWds41zrT8hZktjEQgybmmf2nJ13vGUP9vndnYpCMxS6dSurq6hkQkb+S0BfGNmbU6tGBmLQHdtRUFWtzXlkWPj6Ja2ndsbHgRuzf+FnYkESkkclogWgJfmdn3ZvY9MAtoZ2bfmtmiiKWTHGn9jxTm9v+YmrsWsbZeF9K27Dr+QSIix5HTU0ydI5pCTlrbp7owdfd7nP+fXiytcxl11nxO/CmJxz9QRCQbOR2sT9dSFgDtX7yCyXve5MLB17Gw1uU0XPMpcUnFw44lIgVUTk8xSQGR8vo1TLn6dc7dPJ7UOr3J2KtBd0Ukd1QgCqGOw29mQveXSN7wKfPrXcfBAxlhRxKRAkgFopC6aOQdjO80gBZr32duw5txGQfDjiQiBYwKRCF20fj7GHf+47Rc8Sazm/0X7qALO5KIFCARKxDBwH5fmNkyM1tiZncH68uZ2UQzWxm8lo1UhqLODC7+8mEmNHuQ1gtfZfZ594BTkRCRnIlkCyIduM85VxdoBdxhZvWA/sBk59w5wORgWSLEYoyOc/7NxPr9aD3nBb7q8LCKhIjkSMQKhHNuk3NuQTC/E1gGnA505/dxnIYBPSKVQbyYWKND6nNMOvt2zpv6FLMufSLsSCJSAORLH4SZVQfOBeYAlZ1zm8AXEaBSNsf0NbN5ZjZvy5Yt+RGzUIuNM9oveZkvzuxD67GPMPtPA8KOJCJRLuIFwsxKAh8B/ZxzO3J6nHNukHMu2TmXXLFixcgFLELi4mNos2ww06r0ptVH9zPnhpfDjiQiUSyiBcLMiuGLwzvOuY+D1T+bWZVgexVA41Tno/jEWFqseIuZFbvT8q07mXv74LAjiUiUiuRVTAYMBpY5557LtGkU0CeY7wN8GqkMkrWEUsVosnwEs8t2ptmrtzH/vuFhRxKRKBTJFkQb4Hqgg5mlBlMX4Gmgk5mtBDoFy5LPksoVp96yj1lQqj2Nn7uB1P/3UdiRRCTK5HQ01xPmnJsBWDabUyL1uZJzpSsnctaSUXxbtzP1n7iaRSU+odGDl4YdS0SihO6kLuLKVitJtYVj+C6hMbUeupKlL04KO5KIRAkVCKHCWadQYf541sbXpvrd3Vjx+vSwI4lIFFCBEABOrVeO0rMn8mOx6px+WxdWvTMn7EgiEjIVCDns9HMrEf/lJLbGVqbC9Z1Z+/E3YUcSkRCpQMgRzmx9GhkTJrPTSnNKz078MHZJ2JFEJCQqEPIHZ3U4kz2fTWGfK05i1xQ2Tv0u7EgiEgIVCMlS7S5n8csHk8EdxDqlsHnO2rAjiUg+U4GQbNW/sg4bh00iPn0P+9t2YNvCDWFHEpF8pAIhx9T4+kZ8/+oESh34hZ0tU9i+/KewI4lIPlGBkONq1rcZy58bS4V9P7KtaUd2rt0adiQRyQcqEJIjLe85j0X/Hs1paavZ1LATuzf8GnYkEYkwFQjJsfMeas/ch0Zy5u6l/FD/EtI27ww7kohEkAqEnJAL/n0xX/X7gLN3zGd1nUvZ/+vusCOJSISoQMgJu/B/u/Hlbe9Q99eZLKvTg/Rde8OOJCIRoAIhudJxUC8mXTOEhpsn823tP5GRtj/sSCKSx1QgJNcufucGxl8+kHM3juGbetdwcH962JFEJA+pQMhJueTjvoy9+HmSv/+I+Q364NIzwo4kInlEBUJOWuexd/P5BU/RfOVwvm76Z1zGwbAjiUgeUIGQk2YGl0ztz9jmj9Dy28F83foucC7sWCJyklQgJE+YwcWzHmVcw/tpOfdl5rR7QEVCpIBTgZA8ExNrdFrwDONr3UnL6QOYc8mjYUcSkZOgAiF5KjbOSPn2BSZVv5WW4x9jzhVPhx1JRHJJBULyXFx8DBcsHcgXp11Ly08e5OvrXgg7kojkggqERER8YiytVwxlWqUrafFOP+b1HRR2JBE5QSoQEjEJJeNotnw4M8t1pelrt7Og35thRxKRE6ACIRGVVDaehss+YG7pFBq/cBOL+70ediQRySEVCIm40pUSqLX0U2aXuogGL9zG6r8MCDuSiOSACoTki7Knl6DWsk8ZW7oXZw28nw19HtZ9EiJRTgVC8k3F0+NpsmQ475W+japvPsnmXnfAQQ3LIRKtVCAkX1WpGkubb19lYOkHqPTh/7G9y9WwV8+TEIlGKhCS76qdYVz0zTM8UfpZyox/nz1tOsK2bWHHEpGjqEBIKGrWhJ5f389tpUcQs2Ae+5Nbw6pVYccSkUxUICQ0tWvDXTN60aPUFHb98AsZLVvDrFlhxxKRgAqEhKphQ3hy6nmklJjN+p1lcBdeCG+/HXYsESGCBcLM3jCzzWa2ONO6cmY20cxWBq9lI/X5UnA0bQqvTDibdsVmMTe2NVx/Pdx3H6TrEaYiYYpkC2Io0Pmodf2Byc65c4DJwbIIrVvDm59XICVjAsMr3AXPPQedO6vzWiREESsQzrlpwC9Hre4ODAvmhwE9IvX5UvC0awcfjSrGTTte4LEaQ3AzZkByMixcGHY0kSIpv/sgKjvnNgEEr5Wy2xbr6AgAAA99SURBVNHM+prZPDObt2XLlnwLKOG66CL48EN4fP2N/LnONA7uPwCtWsHrr+vOa5F8FrWd1M65Qc65ZOdccsWKFcOOI/nosstg+HAY/G0LetWcT8Z558Ntt8F118HOnWHHEyky8rtA/GxmVQCC1835/PlSQPTsCUOHwsczK9Oj+DjSH30c3nvPn3JKTQ07nkiRkN8FYhTQJ5jvA3yaz58vBcj118PAgTB6bCy9F/2D9IlTYNcuf8rpP//ROE4iERbJy1zfBWYBtc1sg5ndAjwNdDKzlUCnYFkkW337wvPPw8cfQ5/B7ciYnwopKXDXXXDxxbB+fdgRRQqtuEi9sXPu6mw2pUTqM6VwuvtuSEuDBx+EhISKvDZqNDGvD4J77/V32r38MlxzDZiFHVWkUInaTmqRzPr3h0cegTfegLvuNlzfP/vLX+vX953XvXrB1q1hxxQpVFQgpMB49FH42998g+GBB8CddTZMmwZPPgmffgp16sBbb+lyWJE8ogIhBYYZPPss3HEHDBjgCwaxsf7c04IFcM45cMMNvm9izZqw44oUeCoQUqCYwYsvwi23wGOPwdOHLnNo0ABmzICXXoLZs/3yf/+3xnMSOQkqEFLgxMTAq6/6fukHH4QXXgg2xMb65sXSpdCpkz8Pde65MGVKqHlFCioVCCmQYmNh2DC44gro1w8GDcq0sWpVGDkSPvnE3zeRkgJ/+hN8/31YcUUKJBUIKbDi4uDdd+HSS+H22+HNNzNtNIMePXxr4vHHYexYqFvXd1zs2RNWZJECRQVCCrT4eD+4X0oK3HQTvP/+UTskJsI//gHLl/uC8a9/+c7sQYPUPyFyHCoQUuAlJPgzSm3awLXXwqhRWexUrZpvbkyfDtWrw5//7O+h+PBDXRYrkg0VCCkUkpJg9Gj/dLqePWH8+Gx2PP98f7XTp5/6c1Q9e0KLFjBpkgqFyFFUIKTQKF0axo2DevX82aSpU7PZ0Qy6dYNFi2DIEPj5Z3/VU5s2/g1UKEQAFQgpZMqWhQkToGZN6NoVZs06xs6xsXDjjfDdd/DKK7BhA1xyiW9RjBqlQiFFngqEFDoVK/ozRlWq+Mdaz59/nAMSEuAvf4FVq+C11+CXX6B7d38PxfDhcOBAvuQWiTYqEFIoVani748rW9Y/xvTbb3NwUHw83HorrFjhb7LYt8/3eteoAc88A7/+GvHcItFEBUIKrWrVfJFITISOHf2VrjkSF+fHdFqyBMaM8YMA9u/vb8C7805YuTKiuUWihQqEFGo1a8LkyX4+JQVWrz6Bg2NioEsXf74qNdUPKT5oENSq5ZslH32k009SqKlASKFXu7b/Hb93ry8SP/yQizdp3Nhf8bRunb/ZbvlyP3zHGWfAww9rGA8plFQgpEho2BAmToTt232R2LQpl29UpYp/ctHatf7Gi+bN/ZCyNWv6Ycbffht2787T7CJhUYGQIqNpUz8k06ZNvk9iy5aTeLPYWD8I1KhRvvXwyCP+ctnrr4fKlX0fxsSJkJGRV/FF8p0KhBQprVv7fuc1a/y9cXlyYVK1an4QwNWr/RPurrnGF46LLvLb7r0XvvoKDh7Mgw8TyT8qEFLktGvnR9pYtszfJ7FjRx69cUwMtG3rO7J/+gk++MCfgnrpJX+XdrVq/iqoqVPVspACQQVCiqSLLvLj9C1Y4M8U5Xm3QUKC78T+9FN/Luvtt6FlSxg8GC680Pdl9O3rWxrqs5AoZa4ADCeQnJzs5s2bF3YMKYQ++ACuusr/zv7sM3/PRETt3u07Qj76yHdy79rlb9Br395fUtulix+OXCQPmNl851xyro9XgZCi7q23oE8fPwzTJ5/439f5Yv9+P7Ls55/7adkyv/7ss/25rw4d/PmwcuXyKZAUNioQInlg0CD/iIgrroARI/zN1Plu7Vrfuvj8c/jiC//kOzN/D8aFF/qC0bYtnHJKCOGkIFKBEMkjL7zgn299zTX+8aWxsSGG2b8f5s71hWLKFH8V1L59viO8aVM47zx/SVbr1v5mPbMQw0q0UoEQyUNPPw0PPgg33+wHdo2Jlss49u6F2bN9wZg2Db7++vdna1ep8nuxaN3aj0JbokS4eSUqnGyBCKMhLRK1+veHtDR47DHfYf2f/0TJH+cJCb4ju317v5ye7h94NGvW79PHH/ttMTF+gMGmTX2xaNoUmjSBMmXCSi8FlAqEyFEefdT/cT5ggC8Szz4bJUUis7g4/4u/aVO44w6/7ueffStjwQL45hvf2nj77d+PqVnTF4z69f1j9+rX91dMFS8ezs8gUU8FQuQoZr4opKX5IlGihB+fL+pVruwfdNS9++/rNm/2xeJQ0UhN9ZdqHbqrOzbWF4lDBaNePT+64Vln+We4SpGmAiGSBTN48UV/6v/Q6ab+/cNOlQuVKvlBBC+++Pd1aWl+3KglS2DpUv+6eDGMHHnkcCAVK/pLbrOaypaNwmaV5DUVCJFsxMTAq6/636cPPuiLxN13h50qDyQm+ktnGzc+cv3evb5wrFzpH7+6apUfX2rqVH+zSGalSvmrp7KbTj8dihXLtx9JIkMFQuQYYmP900f37vWXwCYm+hEyCqWEBGjUyE9HS0vz92kcKhzr1vkHa/zwg78cd+vWI/c3g1NPhdNO869VqvjXrOYjfvu65JYKhMhxxMXBu+/C5ZfD7bf736M33BB2qnyWmOj7J+rVy3r7nj2wfr2fDhWOH37wgxZu3Ajz5/v+kKxGtC1d2veflC8PFSr418zzWa1T6yRfqECI5EB8vB8+qWtXuOkmXyR69Qo7VRQpUcJ3bteunf0+GRl+4MKffvLTpk2/v27eDNu2wYYNsHChb5GkpR3780455cSm0qUhKenIKT5efSnHEEqBMLPOwAtALPC6c+7pMHKInIiEBD84a+fOcO21frlbt7BTFSCxsb+fWsqJtDRfNLZu/ePr9u3w22+/T7/+6h/cdGh5796cfUZc3B+LxqGpZMnf50uU8P/gxYv718zTiayLjw/5Fv0Tk+8FwsxigZeBTsAGYK6ZjXLOLc3vLCInKinJP3CoUyfo2dOP1p35AiHJQ4mJULWqn07U/v1HFpDffoOdO/1oujmZfvvNnxo7tLxnjx/q5MCBk/+5zPwpsryc4uL8FBv7+3weDCgWRguiBbDKObcGwMzeA7oDKhBSIJQuDePG+bHzunTRqBbRKR6oGEx5pDjExGdQnH0ksJd4518T2Etxtzf7dewjIZgvxn6KuQN+2n+AuP3BPH6KC16LuQPEsz/T8j7i2HV4v2Iu074cIJYMYl06caQTSwZx+PmTFUaBOB1Yn2l5A9Dy6J3MrC9w6HqRfWa2OB+ynawKwNbj7hU+5cw7FXbtivqMUDC+S1DOvHaMTqHjC6NAZNUj9IcRA51zg4BBAGY272QGnMovypm3CkLOgpARlDOvFaScJ3N8GGNVbgCqZVquCmwMIYeIiBxDGAViLnCOmdUws3jgKmBUCDlEROQY8v0Uk3Mu3czuBMbjL3N9wzm35DiHDYp8sjyhnHmrIOQsCBlBOfNakchZIB4YJCIi+S9anpclIiJRRgVCRESyFHUFwszKmNmHZrbczJaZWWszK2dmE81sZfBaNuSMtc0sNdO0w8z6RVvOIOs9ZrbEzBab2btmlhClOe8OMi4xs37ButBzmtkbZrY58304x8plZg+a2SozW2Fm+XaPdTY5ewbf50EzSz5q/2jK+d/B/++LzOwTMyuTaVu+58wm4+NBvlQzm2Bmp4WZMbucmbb9zcycmVU4qZzOuaiagGHArcF8PFAGeBboH6zrDzwTds5MeWOBn4Azoy0n/qbEtUBisPw+cGMU5mwALAZK4C+cmAScEw05gQuApsDiTOuyzAXUAxYCxYEawGogNsScdfE3Sk0FkjOtj7acFwFxwfwzYX+f2WQsnWn+LmBgNH6Xwfpq+IuA1gEVTiZnVLUgzKw0/oceDOCc2++c244fimNYsNswoEc4CbOUAqx2zq0jOnPGAYlmFof/BbyR6MtZF5jtnNvjnEsHvgQuJwpyOuemAb8ctTq7XN2B95xz+5xza4FV+KFlQsnpnFvmnFuRxe7RlnNC8O8OMBt/b1RoObPJuCPTYhK/39wbVd9l4H+BBzjyBuRc5YyqAgHUBLYAQ8zsGzN73cySgMrOuU0AwWulMEMe5Srg3WA+qnI6534EBgA/AJuA35xzE4iynPjWwwVmVt7MSgBd8H8FRVvOQ7LLldUwMqfnc7aciOacNwNjg/moymlm/zaz9cC1wCPB6mjL2A340Tm38KhNucoZbQUiDt9k+j/n3LnAbnwTPioFN/p1Az4IO0tWgnPj3fFNytOAJDO7LtxUf+ScW4Y/tTARGIdvCp/8SGP5L0fDyESBqMxpZg/j/93fObQqi91Cy+mce9g5Vw2f785gddRkDP64epjfi9cRm7NYd9yc0VYgNgAbnHNzguUP8QXjZzOrAhC8bg4p39EuARY4534OlqMtZ0dgrXNui3PuAPAxcB7RlxPn3GDnXFPn3AX4ZvNKojBnILtcBWUYmajLaWZ9gK7AtS44aU4U5gwMB64M5qMp41n4PwYXmtn3QZYFZnYqucwZVQXCOfcTsN7MDo1AmIIfBnwU0CdY1wf4NIR4Wbma308vQfTl/AFoZWYlzMzw3+cyoi8nZlYpeD0DuAL/vUZdzkB2uUYBV5lZcTOrge9o/zqEfMcTVTnNP0Ds70A359yeTJuiJqeZnZNpsRuwPJiPmozOuW+dc5Wcc9Wdc9XxRaFp8Hs1dznzo7f9BHvmmwDzgEXASKAsUB6YjP+rcjJQLgpylgC2AadkWheNOf+F/495MfAW/iqGaMw5Hf/HwEIgJVq+T3yh2gQcCP6Hu+VYufBN/NXACuCSkHNeHszvA34GxkdpzlX48+OpwTQwzJzZZPwo+H9oEfAZcHo0fpdHbf+e4Cqm3ObUUBsiIpKlqDrFJCIi0UMFQkREsqQCISIiWVKBEBGRLKlAiIhIllQgRHLJzB4ORks9NMpny7AzieSlfH/kqEhhYGat8Xf+NnXO7QuGVY4POZZInlKBEMmdKsBW59w+AOfc1pDziOQ53SgnkgtmVhKYgb+jfhIwwjn3ZbipRPKW+iBEcsE5twtoBvTFD1E/wsxuDDWUSB5TC0IkD5jZn4A+zrnLws4iklfUghDJBfPPJc88wmcT/CMeRQoNdVKL5E5J4D9mVgb/kJtV+NNNIoWGTjGJiEiWdIpJRESypAIhIiJZUoEQEZEsqUCIiEiWVCBERCRLKhAiIpIlFQgREcnS/wfT5idyp1X72gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "BS.plot([60,140,0,50])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec1'></a>\n",
    "## Binomial tree\n",
    "\n",
    "One of the most used methods for pricing American options is the Binomial tree. \n",
    "\n",
    "We have already encountered the binomial tree in the notebook **1.1**. The following algorithm is almost a copy/paste from that notebook.     \n",
    "There are just two additional lines:\n",
    "- `S_T = S_T * u`. This an efficient method to retrive the price vector at each time steps.  \n",
    "- `V = np.maximum( V, K-S_T )`. This line computes the maximum between the conditional expectation V and the intrisic value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "S0=100.0    # spot stock price\n",
    "K=100.0     # strike\n",
    "T=1.0       # maturity \n",
    "r=0.1       # risk free rate \n",
    "sig=0.2     # diffusion coefficient or volatility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "American BS Tree Price:  4.81624866310944\n"
     ]
    }
   ],
   "source": [
    "N = 25000              # number of periods or number of time steps  \n",
    "payoff = \"put\"        # payoff \n",
    "\n",
    "dT = float(T) / N                             # Delta t\n",
    "u = np.exp(sig * np.sqrt(dT))                 # up factor\n",
    "d = 1.0 / u                                   # down factor \n",
    "\n",
    "V = np.zeros(N+1)                             # initialize the price vector\n",
    "S_T = np.array( [(S0 * u**j * d**(N - j)) for j in range(N + 1)] )  # price S_T at time T\n",
    "\n",
    "a = np.exp(r * dT)    # risk free compound return\n",
    "p = (a - d)/ (u - d)  # risk neutral up probability\n",
    "q = 1.0 - p           # risk neutral down probability   \n",
    "\n",
    "if payoff ==\"call\":\n",
    "    V[:] = np.maximum(S_T-K, 0.0)\n",
    "elif payoff ==\"put\":\n",
    "    V[:] = np.maximum(K-S_T, 0.0)\n",
    "\n",
    "for i in range(N-1, -1, -1):\n",
    "    V[:-1] = np.exp(-r*dT) * (p * V[1:] + q * V[:-1])    # the price vector is overwritten at each step\n",
    "    S_T = S_T * u                    # it is a tricky way to obtain the price at the previous time step\n",
    "    if payoff==\"call\":\n",
    "        V = np.maximum( V, S_T-K )\n",
    "    elif payoff==\"put\":\n",
    "        V = np.maximum( V, K-S_T )\n",
    "    \n",
    "print(\"American BS Tree Price: \", V[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec2'></a>\n",
    "## Longstaff - Schwartz Method\n",
    "\n",
    "This is a Monte Carlo algorithm proposed by Longstaff and Schwartz in the paper [1]:\n",
    "[LS Method](https://people.math.ethz.ch/~hjfurrer/teaching/LongstaffSchwartzAmericanOptionsLeastSquareMonteCarlo.pdf)\n",
    "\n",
    "The algorithm is not difficult to implement, but it can be difficult to understand.  I think this is the reason the authors started the paper with an example.    \n",
    "Well. I think they had a good idea, and for this reason I want to reproduce their example here. \n",
    "\n",
    "The same code is copied in the class `BS_pricer` where the function `LSM` is implemented.\n",
    "\n",
    "**If in the following code you feel that something is unclear, I suggest you to follow the steps (not in python) proposed in the original paper**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 4          # number of time steps\n",
    "r = 0.06       # interest rate\n",
    "K = 1.1        # strike \n",
    "T = 3          # Maturity\n",
    "\n",
    "dt = T/(N-1)          # time interval\n",
    "df = np.exp(-r * dt)  # discount factor per time interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = np.array([\n",
    "            [1.00, 1.09, 1.08, 1.34],\n",
    "            [1.00, 1.16, 1.26, 1.54],\n",
    "            [1.00, 1.22, 1.07, 1.03],\n",
    "            [1.00, 0.93, 0.97, 0.92],\n",
    "            [1.00, 1.11, 1.56, 1.52],\n",
    "            [1.00, 0.76, 0.77, 0.90],\n",
    "            [1.00, 0.92, 0.84, 1.01],\n",
    "            [1.00, 0.88, 1.22, 1.34]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1.0 & 1.09 & 1.08 & 1.34\\\\1.0 & 1.16 & 1.26 & 1.54\\\\1.0 & 1.22 & 1.07 & 1.03\\\\1.0 & 0.93 & 0.97 & 0.92\\\\1.0 & 1.11 & 1.56 & 1.52\\\\1.0 & 0.76 & 0.77 & 0.9\\\\1.0 & 0.92 & 0.84 & 1.01\\\\1.0 & 0.88 & 1.22 & 1.34\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1.0  1.09  1.08  1.34⎤\n",
       "⎢                     ⎥\n",
       "⎢1.0  1.16  1.26  1.54⎥\n",
       "⎢                     ⎥\n",
       "⎢1.0  1.22  1.07  1.03⎥\n",
       "⎢                     ⎥\n",
       "⎢1.0  0.93  0.97  0.92⎥\n",
       "⎢                     ⎥\n",
       "⎢1.0  1.11  1.56  1.52⎥\n",
       "⎢                     ⎥\n",
       "⎢1.0  0.76  0.77  0.9 ⎥\n",
       "⎢                     ⎥\n",
       "⎢1.0  0.92  0.84  1.01⎥\n",
       "⎢                     ⎥\n",
       "⎣1.0  0.88  1.22  1.34⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_matrix(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous cell we defined the stock matrix S.   \n",
    "It has 8 rows that correspond to the number of paths.  \n",
    "The 4 columns correspond to the 4 time steps, i.e. each row is a path with 4 time steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Example price=  0.11443433004505696\n"
     ]
    }
   ],
   "source": [
    "H = np.maximum(K - S, 0)           # intrinsic values for put option\n",
    "V = np.zeros_like(H)               # value matrix\n",
    "V[:,-1] = H[:,-1]\n",
    "\n",
    "# Valuation by LS Method\n",
    "for t in range(N-2, 0, -1):\n",
    "\n",
    "    good_paths = H[:,t] > 0        # paths where the intrinsic value is positive \n",
    "                                   # the regression is performed only on these paths \n",
    "    \n",
    "    rg = np.polyfit( S[good_paths, t], V[good_paths, t+1] * df, 2)    # polynomial regression\n",
    "    C = np.polyval( rg, S[good_paths,t] )                             # evaluation of regression  \n",
    "    \n",
    "    exercise = np.zeros( len(good_paths), dtype=bool)    # initialize\n",
    "    exercise[good_paths] = H[good_paths,t] > C           # paths where it is optimal to exercise\n",
    "    \n",
    "    V[exercise,t] = H[exercise,t]                        # set V equal to H where it is optimal to exercise \n",
    "    V[exercise,t+1:] = 0                                 # set future cash flows, for that path, equal to zero  \n",
    "    discount_path = (V[:,t] == 0)                        # paths where we didn't exercise \n",
    "    V[discount_path,t] = V[discount_path,t+1] * df       # set V[t] in continuation region\n",
    "    \n",
    "V0 = np.mean(V[:,1]) * df  # discounted expectation of V[t=1]\n",
    "print(\"Example price= \", V0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix `H = np.maximum(K - S, 0)`, is the matrix of intrinsic values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.1 & 0.01 & 0.02 & 0.0\\\\0.1 & 0.0 & 0.0 & 0.0\\\\0.1 & 0.0 & 0.03 & 0.07\\\\0.1 & 0.17 & 0.13 & 0.18\\\\0.1 & 0.0 & 0.0 & 0.0\\\\0.1 & 0.34 & 0.33 & 0.2\\\\0.1 & 0.18 & 0.26 & 0.09\\\\0.1 & 0.22 & 0.0 & 0.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0.1  0.01  0.02  0.0 ⎤\n",
       "⎢                     ⎥\n",
       "⎢0.1  0.0   0.0   0.0 ⎥\n",
       "⎢                     ⎥\n",
       "⎢0.1  0.0   0.03  0.07⎥\n",
       "⎢                     ⎥\n",
       "⎢0.1  0.17  0.13  0.18⎥\n",
       "⎢                     ⎥\n",
       "⎢0.1  0.0   0.0   0.0 ⎥\n",
       "⎢                     ⎥\n",
       "⎢0.1  0.34  0.33  0.2 ⎥\n",
       "⎢                     ⎥\n",
       "⎢0.1  0.18  0.26  0.09⎥\n",
       "⎢                     ⎥\n",
       "⎣0.1  0.22  0.0   0.0 ⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_matrix(H.round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix V contains the cash flows.\n",
    "\n",
    "**Important**    \n",
    "To simplify the computations, the discounted cashflows are reported at every time steps. \n",
    "\n",
    "For instance:    \n",
    "In the third row the final cashflow (0.07) is discounted at every time step, till t=1.    \n",
    "In the paper, the authors just consider the cashflow (0.07) at time t=3 and the discount is performed at the end of the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0621 & 0.0659 & 0.07\\\\0.0 & 0.17 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.34 & 0.0 & 0.0\\\\0.0 & 0.18 & 0.0 & 0.0\\\\0.0 & 0.22 & 0.0 & 0.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0.0   0.0     0.0    0.0 ⎤\n",
       "⎢                         ⎥\n",
       "⎢0.0   0.0     0.0    0.0 ⎥\n",
       "⎢                         ⎥\n",
       "⎢0.0  0.0621  0.0659  0.07⎥\n",
       "⎢                         ⎥\n",
       "⎢0.0   0.17    0.0    0.0 ⎥\n",
       "⎢                         ⎥\n",
       "⎢0.0   0.0     0.0    0.0 ⎥\n",
       "⎢                         ⎥\n",
       "⎢0.0   0.34    0.0    0.0 ⎥\n",
       "⎢                         ⎥\n",
       "⎢0.0   0.18    0.0    0.0 ⎥\n",
       "⎢                         ⎥\n",
       "⎣0.0   0.22    0.0    0.0 ⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_matrix(V.round(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] F. Longstaff, E. Schwartz (2001) \"Valuing American Options by Simulation: A Simple Least-Squares Approach\", The Review of Financial Studies, vol 14-1, pag 113-147.   "
   ]
  }
 ],
 "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
}