LSMC2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate some price paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 1\n",
    "N = 3\n",
    "dt = T/N\n",
    "M = 200\n",
    "S = np.zeros(shape=(M, N+1))\n",
    "S[:, 0] = 1  # random starting price\n",
    "\n",
    "k = 1.1\n",
    "\n",
    "mu = 0.06\n",
    "si = 0.3\n",
    "\n",
    "# for i in range(1, N+1):\n",
    "#     dz = np.random.normal(size=M) * np.sqrt(dt)\n",
    "#     dS = ((mu * dt) + (si * dz)) * S[:, i-1]\n",
    "#     S[:, i] = S[:, i-1] + dS\n",
    "\n",
    "# df = pd.DataFrame(S)\n",
    "\n",
    "# df.transpose().plot(color=\"red\", alpha=0.3)\n",
    "# plt.legend([])\n",
    "# plt.plot([0,N], [k, k], label=\"strike price\")\n",
    "# plt.show()\n",
    "\n",
    "# df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABBg0lEQVR4nO29e3BkV37f9z0AGsAAg8fgjcFjMO8XhsOZATnkklLoyKpQ8lprRVuVXSVyVmXXuJSsHamSKjuPEmWlUrLKTipKNjGL2TAUXaldS5YsU1uryK5kKe7WLrkczvI5HM6LGLwHbzTQ6Pc9+eOLs7e70S8At/s++vep6gLQfdF9bt97v/d3fuf3UFprCIIgCP6nzu0BCIIgCM4ggi4IghAQRNAFQRACggi6IAhCQBBBFwRBCAgi6IIgCAGhpKArpV5VSi0qpT4uss0LSqn3lVKfKKX+ytkhCoIgCOWgSsWhK6V+FsAWgNe11uN5Xu8E8EMAL2qtp5RSfVrrxUoMVhAEQShMSQtda/0WgNUim/wqgD/VWk/tbC9iLgiC4AINDrzHGQAhpdSbANoA/IHW+vV8GyqlbgC4AQCtra3Xzp0758DHC4Ig1A7vvffesta6N99rTgh6A4BrAH4OwCEAP1JKva21vpu7odb6FQCvAMDExIS+efOmAx8vCIJQOyilHhV6zQlBnwGworWOAIgopd4CcBnALkEXBEEQKocTYYv/BsDzSqkGpVQLgOsAPnXgfQVBEIQ9UNJCV0p9C8ALAHqUUjMAXgIQAgCt9cta60+VUv8PgA8BWAC+qbUuGOIoCIIgVIaSgq61/moZ2/xTAP/UkREJgiDUCMlkEjMzM4jFYrtea25uxvDwMEKhUNnv54QPXRAEQdgHMzMzaGtrw9jYGJRSP31ea42VlRXMzMzg+PHjZb+fpP4LgiC4RCwWQ3d3d5aYA4BSCt3d3Xkt92KIoAuCILhIrpiXer4Y4nIRhKASiwGzs0BDA9DUBDQ389HYCNSJLRdERNAFIYgkEsAPfwhEIvlfb2ykuGcKvfk982d9fXXHLRwIEXRBCBqpFPDOO0A0CnzhC0BLCxCP02LP93Nzkz/zFeoLhbIFvtANYA+RGEI2Wuu87pVShRPzIYIuCEHCsoCbN4GNDeCpp4Dubj5/6FDx/9OaVn0x4V9b4++Wtfv/6+qKW/rmZ2MjsA/fcFBpbm7GysrKroVRE+XS3Ny8p/cTQReEoKA18JOfAEtLwJNPAv395f+vUrYF3t5efNtksrTFv7zM7Yp9TjGXT1NTTfj5h4eHMTMzg6WlpV2vmTj0vSCCLghB4ZNPgLk54Px5YGSkcp8TCvFx+HDx7dLp4sIfjQLr6/w9H7l+/kI//eTntyzgs894s+3qQigU2lOceSlE0AUhCNy9C3z+OXDyJHDqlNujIfX19N+3tBTfzrLo7ikk/LFYcT9/bhRPIeF328+/scEZ1OYmv5uuLsc/QgRdEPzOo0e0+oaHaZ37jUz/ezG0phunmPCvr/NnOl34c4ot7pqwTif9/FoD9+/zGDU1AdevA319zr1/BiLoguBn5ueBDz/kFP7y5WAvOCpFsW1sLL1tKlVc+CMRYGVlb37+fD9L+fm3tmiVr68DQ0PApUsVnSmIoAuCX1leBm7d4tT92rWaWEQsm4YG+vhL+fktq7jwG6s/kSgc1lloUXdhAZicZITRtWvA0aOV2NMsRNAFwY9sbADvvgu0tgJPP+2vhUEvUVdXnp9fa4p8sUXerS17sffhQ/rKOzuB48eBDz4A7tyxBX9wsCICL4IuCH4jEmHiUCgEPPOM+4t9tYBStiXe0VF4u6kp4P33uTh98iTQ07Nb+Dc2ir/HARBBFwQ/EYsBb79Ni/GZZ0ovJArVIRajFb64CPT2Mg+gVDJXBSjpdFNKvaqUWlRK5e1CpJR6QSm1oZR6f+fx284PUxAEJJO0zONxRkqU8g8L1WF2FnjzTS6yjo/zRuuCmAPlWeivAfgGgNeLbPN9rfUXHRmRIAi7SafpM9/aos+8s9PtEQmJBCOM5ueBI0eAK1e4puEi5bSge0spNVaFsQiCkA+tGc2ysgJcvcopveAuCwsU82SSsf8nT3oiZNQpH/qzSqkPAMwB+K+01p/k20gpdQPADQAYHR116KMFIeB8+CEFZHycscyCeySTLLEwPc2aN888U7r2TRVxQtBvATimtd5SSv0igD8DcDrfhlrrVwC8AgATExN7rw0pCLXGnTuMnDh9muFvgnssLzOCJRbj8ThzxnOx/wcejdY6rLXe2vn9uwBCSqmeA49MEGqdhw+Be/eAY8eAc+fcHk3tkk4DH30E/OhHjPd//nkeD4+JOeCAha6UGgDwWGutlVJPgzeJlQOPTBBqmdlZTu0HB5kuLrjD2hpT9yMR4MQJCrmHk7hKCrpS6lsAXgDQo5SaAfASgBAAaK1fBvBlAL+hlEoBiAL4it5Pqw1BEMjiIkWku5uLoB5YbKs5TJnbBw8Y6/+FL9jNQjxMOVEuXy3x+jfAsEZBEA7K2ho7DrW3MzzRg9P6wBMOM6pocxMYHQUuXmRtGB/gj1EKQi2wucnEoeZmJg75REQCQxXL3FYKOWMEwQtEo0zpr6tjKFxTk9sjqi2qXOa2UoigC4LbJBIU81QKeO650pX/BOfQmp2ePv2UM6IqlbmtFCLoguAm6TTw4x8D29ueS1IJPNvbjCtfWbEbhPh8ZiSCLghuYVmsz7K+DkxM+CKKIjBMTTEsFGBlxEo21a4iIuiC4AZa0zpcWqJlODDg9ohqg8wytz09rpW5rRQi6ILgBrdvM3no3DmGxgmVZ3aWGZ+Wxbo4Y2OBi/EXQReEanPvHtP6T5xgTRChsiQSFPK5Oc+Uua0UIuiCUE2mplhwa2gIuHDB7dEEn8eP6WLxWJnbSiGCLgjVwtTQ7uuj7zbAwuI6Hi9zWylE0AWhGqysAO+9x05DExOS0l9JfFDmtlKIoAtCpQmHGWve0sL6LB6u1udr0mkmCH3+OX3kzz1Hn3kNIYIuCJVke5tZoA0NnPY3Nro9omDiszK3lUIEXRAqRTxOMbcsWosBinf2DLllbp99lvHlNYoIuiBUgmSSYh6LUWTa2tweUfAIh2mVh8O+K3NbKWp77wWhEpiU/s1N+sxrzI9bcUyZ27t3WRHx6adZi0Uo3VNUKfWqUmpRKfVxie2eUkqllFJfdm54guAztGZzhJUVJrD4rJ6259naAn7wA8byDw4CL7wgYp5BObE8rwF4sdgGSql6AL8P4N86MCZB8C8ffQTMz3P6PzTk9miCg9bMrv2rv+LC57VrbM8ni8xZlNOC7i2l1FiJzf4+gD8B8JQTgxIEX3LnDvDoEWOfT5xwezTBIYBlbivFgX3oSqkhAL8M4K+hhKArpW4AuAEAo1KQSAgSn3/OGi2jowyZE5whs8zt5ctSyKwETiyK/s8A/qHW2lIlUpm11q8AeAUAJiYmtAOfLQjuMzsLfPwxS+A+8YTbowkGAS9zWymcEPQJAN/eEfMeAL+olEpprf/MgfcWBG+ztER3QHc3fbpSn+Xg1ECZ20pxYEHXWh83vyulXgPwHRFzoSZYX2d44uHDwFNP1WRmoqPUUJnbSlFS0JVS3wLwAoAepdQMgJcAhABAa/1yRUcnCF5lawt45x0uzl2/7ssO8Z4is8ztuXPAqVNile+DcqJcvlrum2mtv3ag0QiCH4jFmAWqFOuzNDe7PSL/UqNlbiuFZIoKwl4wKf3JJPCFL4hL4CDUcJnbSiGCLgjlkk7TzRKJ0M3S0eH2iPyJlLmtGCLoglAOlgXcvMkyrRMTNV3R70Bklrk9fpxt4WQx2TFE0AWhFFrbMdFPPMEaIsLekDK3VUEEXRBKcfs2MDMDnD0LHDvm9mj8h5S5rRryrQpCMe7fZ1Go48e5aCeUj5S5rToi6IJQiOlpLt4NDdGqFMpna4tW+fo6cPQocOmSVEasAiLogpCPhQX6zXt7WUdEklzKQ2tGr3z6KRc7r12joAs2Kyu8uVWgi5UIuiDksroKvPcewxInJiQ2ulxyy9w+8YQkXQG8ya2vs07+T37CxeGJCeCXfsnxjxJBF4RMwmHgxz9mZb/r12XxrlykzG02ySSjoh4/ZgG31VVgcpLGwZUrTEqrAHK2CoJhe5uJQ/X1DKsTn29pYjHgww8pXN3ddE+1tLg9KncIh20RX1ujZV5fz5j7eJzfzZUrLLNcIUTQBQHgBff228xifO45qb1dDqbMbTpdm2Vu02mWL3j8mEIejfL5jg4WFwM4c2luZgLVuXMVn/GJoAtCKkXLPBZjcagKLFYFitwyt08+yRLCtcD2ti3gy8tMmGpoYJLUmTN2U/CPP6bPvL2d4ZqdnVUZngi6UNtYFmuah8Osad7V5faIvE2tlbm1LPq/jStla4vPt7ZyRtLXR1dTXR1dLI8eMcLHsmiVnzhR1UV1EXShdtEauHWLltaVK5L0UoxUioueU1PBL3Mbj2cvaKZSFOXubmYK9/fvrrK5ucm1hNVVhrpeuuRKJc5yGly8CuCLABa11uN5Xv8SgP8egAUgBeA3tdY/cHqgguA4Zlp84QIwPOz2aLxL0Mvcag1sbNiulPV1Pt/czBj6/n66VPL5vy2LzcHv3+frV664ei6VY6G/BuAbAF4v8Pr/C+ANrbVWSj0B4I8ASNtzwdvcvcswspMn+RB2E+Qyt8kkre/FRT7icT5/5AhdSf39pWcgKyt0P0UiFPGLF12PjCqnY9FbSqmxIq9vZfzZCkA7MC5BqByTk0zuGBmhdS7sJohlbjc3bVfK6iot81CIfvD+frpKyhHkZJIF26amGKL5zDP8Xw/giA9dKfXLAH4PQB+Av1FkuxsAbgDAaK0nHgjuMDfHCI3+fibACNlYFmcv9+/7v8xtOk0r2rhStrf5fHs7Z2X9/bTI97KoOzdHV10iwQXhM2c8daNzRNC11v8awL9WSv0s6E//6wW2ewXAKwAwMTEhlrxQXZaWaHV2dbHGSJCjM/ZDbpnbCxf81/w6Gs0OK0ynKbg9PRTgvr795RhEo1z0XFxkCKJHF4UdjXLZcc+cUEr1aK2XnXxvQTgQ6+vsONTayrhgD1lVruPnMrdaZ4cVbm7y+ZYW3pT6+ijm+13ENcXG7tzh3xcv0gXlUWPgwIKulDoF4MHOouhVAE0AVg48MkFwikiEiUOhEC0rv1mdlWRrixEsa2v+KXObSGSHFSaTFNjubq6L9Pc7k+i0scFFz40NvuelS57PIC4nbPFbAF4A0KOUmgHwEoAQAGitXwbwKwD+tlIqCSAK4D/SWos7RfAGsRjwox/x92eflep/Br+Vuc0MK1xb43NNTayLYhY0nUqrT6e5aP7wIW9uXv9uMignyuWrJV7/fQC/79iIBMEpkknWZ0kmKeYuJHp4kswyt319XBz22o0ulcqukxKL8fnOTrYC7OtjzRSnXR9LS/SVb28ziej8eV/N6CRTVAgm6TTL4EYiVa2l4Xm8XOY2EqGAm7BCUyelr89+NDVV5rPjcX4vs7N01zz3nC/LQIigC8FDazaoWF3ldNkjMcKu4sUyt5bFWYLxh0cifL6tjQuPJqyw0lmp09MU83Sa1v+pU77NhBVBF4LHBx9QIC5d8o3vs6LMzVHM02n3ozRisewFzXSa4tnTw0JWfX3Vu9FEIvxelpdpjV++7PuqkSLoQrD49FNaXGfOsBperWLqkzx4QEHv7GSdkWoLltZcxDQiHg7z+UOHGJFiwgqrGUZqWfxe7t7l5z7xBF1PHg1F3Asi6EJwePCA8dTHjnHqXEtYFmPtV1boalpdtasEVrvMbW77tUSCn93VxUXG/n73as6vrXEGt7nJ2dvFi95bED4AIuhCMJiZYX2NwUG6WoJOOk1xWl2liK+t8TmAYjk8TF95d3flFhIzCYezwwq1ZshfZp0UN6NFUinO3iYnOTvwU/LUHhBBF/zP48cMw+vpAa5eDcTUeReplG15r6zQGrcs7mt7O2cl3d20gquRGJRO29UKHz+2wwo7Olhit6+Pbh4vHIuFBdbvicXopz97NrDNv4O5V0LtsLrKiJb2dnYc8ml0wi6SSdt9srJCf7jWFMjOTgpTdzejQKpl+Zr2a48fc0wmrLC31w4r9JL7IhajkC8s2OdHwMNXRdAF/7K5yVjz5mbg+nV/W13xeLaAm8XDujqK9unTtoBXawHRtF8zrpTc9mv9/ZwReO0m6oFWcG7h4ytAqGmiUWaB1tWxPks1/MROEo1SLJeX+dOIZX09RfLcOQp4Z2d1hahY+zXTQ9PLGbceaQXnFiLogv9IJFifJZ0GvvAF9xNkyiESybbATW3uUIgCPjpK0Wxvr66AZ7Zfe/yYvwPltV/zEqaO+4MHnmgF5xYeP0qCkEMqxcqJ0Sjrs3iwJjUAWoqZAm4WDRsbKdzGB97WVv2Fw3zt15TaW/s1L+HBVnBuIYIu+AfLYk3zjQ1gYsI7tTa0ps97ZcUW8USCrzU32+GD3d3uZSI61X7NS3i4FZxbiKALB8eygPl5+n/7+irjMtCaoYlLS6xDMjDg/GeUi2XxppKZxJNM8rWWFgqkCSF0038bjTLRKrf9muncs9f2a15idpb1VzzaCs4tRNCFgzE/z2gCU1gpFAKGhjj1dbJDvKmEd/48U8ariWVlJ/GsrtpJPIcP09dsBNwrDRBWV4F33+U4e3spev393gor3A/b2wxF9HgrOLcQQRf2x/IyhXx9nX5gEwM+M8NaKpOTtE6Hh/k4yMLlvXtsxnDiBIWp0qTT2eK9tkZRBygeZgGzq8ub0TVTUxS9Q4dYBtbnBacAZLeCUwoYH2fUjV9nGBWinI5FrwL4IoBFrfV4ntf/YwD/EIACsAngN7TWHzg9UMEjbGxQyJeWKBhPPknBNhdWXx8XLufnKe6ffcZHVxct68HBvSXCPHrEi3h4mE2LK0EymS3g6+t2Ek9HB6sTGgH3crMDrelTfviQVvm1a94eb7n4sBWcW5Rjob8G4BsAXi/w+ucA/j2t9ZpS6hcAvALgujPDEzxDJEJhnp2lSFy8SAspn7+8oYHiPTJCP+7MDB8ffEDLcWCAAt3bW9zfPj/P7U1XHaessUTCFu/l5ewkns5OzgJMEo/Xw/UMySRw6xZdEceP8/j43XrNbQU3MUGDQChIOS3o3lJKjRV5/YcZf74NoPaCP4NMPM743kePKHinTwMnT5Zv+R06xP85fZqW78wMbwpzc7xIjb89NyV7eZkC1dnJC/kgC62xWHYEiukMX19P0T571k7i8ePCWiRid2d64gnWdfE7i4tMEIpGfdkKzi2cNj/+DoC/KPSiUuoGgBsAMOql1lfCblIpJmk8eED/8egoIwkOsqjW2cnHhQt02UxP80bx+ef08xp/eyLBBb3WVqb071Vkt7ezBdws2DY00G1iKhF2dPg/HXx5maGcSjEuv7vb7REdjIC0gnMLxwRdKfXXQEF/vtA2WutXQJcMJiYmtFOfLTiIZXFB8949CuvRo0w2cTL8rq6OvtD+froK5uZoud+5w9DE2Vm6WX7mZ8qzyra2sgU8GuXzoZCdsm6yMP3uhshkchL4+GMK39NP+yNjthhTU1wDCEArOLdwRNCVUk8A+CaAX9BarzjxnkKV0ZpC+tlntHB7eynkla5OFwpxSn3sGMX4z/+cVlpDA/DWW/S3j4ww/VwpjtNkYRoBj8f5Xk1NFG7jAz98OFgCbrAsWrGTk7wpXr3qH19/PiIRrq+srASmFZxbHPgsUEqNAvhTAL+mtb578CEJVefxY1rH4TDdEG5k3CWTXADt7QW+9CUK9/S0fZOJx+mPb2iwMxoPHeL2JguzFoowJRIsF7y8zBvXuXP+vWnltoK7fJk3b7/ujwcoJ2zxWwBeANCjlJoB8BKAEABorV8G8NsAugH874oHIqW1nqjUgAUHWVtjCOLKCsXw2jVGEVT7gkqn6TPf2uICqOkGv73N3xMJWuKxGC23wUH64U+e9H+izF7Y2uLiZzTq/+JTua3gxse9GdPvM5TW7riyJyYm9M2bN1357Jpna4tCvrDAi+jMGS56uuGvTKWA732PKepDQ7TUTBJPW5ttfZu628bfvrbGG09PD626gQF/RqiUy+IiLfP6eiZxOZmFW02SSc4GTSu4S5cC2Qqukiil3itkNPvY8SbsmWiU09vpaQrDuXOMWa6m/9W0UjP+bxM7fewY/fWZAp6vWNTYGB+RiO2SuXWL+zA4aEewBGna/uABb8BtbVz89GtSjWkFF48HvhWcW8i3WQskk3b6PEARP326OtX1jLvELGKGw3YW5sYGbywvvkih2kuccWsrb0hnz/L9Z2ZovU9P0w1jQiDd6i7vBJbFWOzpad6srlzx5yykBlvBuYUIepBJpyni9+9T1EdGKICVtPBMKzVjgRdqpba2xqn3+DgXw/aLUrZVPz7OBd6ZGVq19+9zkXd4mO4cP/lo43HGl6+u0iV25oz/Zh2ZreC05rrH8eMSilhBRNCDiNaM6b17l9ZRfz+t2UpVpdvYoBW5uGgn8RRrpTY7SzEfHGRmo1PU13OB7ehRCuLsLMX9k08Y39zXR3Hv7/e2pRsOc/EzHudC9dGjbo9o72xuctFzba0mW8G5hQh60MgsZ9vVRUGoRKZdMknBnJqioNfV8cI9dqx4K7XFReAnP+E2V69WzupsaqKf9sQJioupJ/P4Mf22R49S3Lu6vGX5LixwTSAUYpak31wTphXc/fvcB79H4/gMEfSgkK+cbSWaQKysUMTn5njxtrfT+hoaKu0DX1ujGyGz3G41aGtjLZBz5zh+U0/GdLoZGqI7ym0L8t49zlw6O/n9+C0kc3mZPv9IhN/nhQv+64Lkc0TQ/U6pcrZOEI/TpTI1xYu1oYFhjqOj9FGXw9YWe4E2NTFxyY1CSybMsaeHN6GFBe7X/fsU0yNH+N0dPVpdIUqn6Z6YneXN5fJlb7uEckkk6NKanuZN8dln+R0LVUcE3a/klrN1esFJa7pHpqboptCabpIzZ+j73ovgRKPAj37EsT3zjDcWJ+vrKZ5DQ1xnMP72jz6izz3T317JmUQsxqSq9XXOIE6frtxnVYLZWdaTSSbtqpp+uhkFDBF0vxGP05qcnNxfOdtSbG9TxKenKTZNTXz/kZH91ddIJIC332b8+XPPue/WyEdzM/fx5EkuSBp/+8ICv9ejR7n/TifzrK9TzJPJyrnIKkVuK7jLl6UVnAcQQfcLlShnazBNnqem6AdVihbqpUsHa/qcTjNaY3vbP70f29s52zl/nt/F9DTF/dEj51rqAbRs33+fN8znn/fHdwNwpvbwIWeH0grOc4ige51KlrMNhyniMzO0Elta+N4jIwe/UVgWF0DX11mfxW91upVi1E5vb+GWesbfvpfZkdb8/3v3+B4TE95wQZVDZiu4gQGKuV+zVgOKCLpXyS1n29NDq/GgYWyplB3hsb5O63tw0G587ISlpTWtz8VFTsX95ErIR6GWeh9+SP9xfz9fK9VSL51myOb8PL/vS5f8kWSTTjP65vPPefORVnCeRQTdiywuMnLFyXK2q6t2uGE6zVC+8XEuCjod0XH7Nm8a585RuIJEoZZ68/PFW+pFo3Q/bW6y3+eJE26Mfu9IKzhfIYLuJXLL2V69yin9fq3meJyCMzXFsMGGBgqOKYRVCe7fp4/V1IsJMrkt9YyvPbelXjRK91M6zZo1fX1uj7w00grOl4ige4HccraXLu2/nK3WFBcTbmhZvBCffJLT5EpWt5ua4n4MDdEKrRVyW+rNz3Mx9c4d4Pvf5+xoeJhFyPyQ+Smt4HyLCLqbxGL0kTtRzjYatcMNo1FO/8fGeGOodMVBy6Ir58MP6Rp68snajXoIhfidj4wwhX92ljfppibghz/kesLwML8nr31Hma3gurtZZ0dawfmKcjoWvQrgiwAWtdbjeV4/B+D/AnAVwH+rtf5njo8yaDhVztayaNVPTdEqBygUFy9WPiEmHGZY39ISBSCdpvU5MSHWXCrFZhSLi3RVXLzIyBDjbzcib0oOuB2yKK3gAkM5puBrAL4B4PUCr68C+AcA/pYzQwowueVsh4dple819Gtz0w43TCT4/2fP8iKsVBhZLEbxXlqikJvGzK2tdoRHb69kCUYidju9S5c4SwKYlHTkCMV9cZEzqclJrje0t9slfqtdv2V1lTMraQUXCMpqQaeUGgPwnXwWesY2vwNgq1wLfb8t6P7xn3+C23PhPf+f68Tjdo/MxkbGfO9F/LSmeMditACV4vS+ubkyUQdaA6kkkEjy5pNO8/k6BYQa+ZmhkFjjmaSSQHiTv7e1lT4uWvO8iMd5TAH+T1MT0NQIoIIWstY8H2MxHsPDrTyuQlW4cLQdL/3N/a0zeaYFnVLqBoAbADAatHC2QiQSvHDSafrG29r25iNPpYB4DIgneBHW19MqbmpyfkqcSlG8kwkguSMwSnG8zU284GvdAi9EPAZsRfj9tLcBdWV8T0rxhtzcDFhpW9y3toCI4o2/qcn5G3YiwZmEZQGHmoFDLeJeCQhVFXSt9SsAXgFooe/nPfZ7V6s6Py1nG6aInztXfoJNImGHG25u7hSSGqZrw8nQsa0t24WyvGxbiZ2ddJ/09NjNmYsRjzPkMpXituahVPbf+R6Z2/gRrRne9/kywxGvXj2YAGvN73J6mtEyyQ3nWur9tBXcMt08ly/7I+pGKBuJcnGacJghX3stZ6s1RXVqigudlmUXPTp61Jlww3jcXshcWuIFDtg1wXt7Gd1QbHHWsrjAt7ZmP6LRg48NKC72pW4G+93mIDedZJKLn0tLLOx1/vzBLV2leBPt6nKupZ60gqsZRNCdYnubccd7LWcbi9nhhtvb/N+xMWeiH9JpRqAYK9z09wyFaH2bhcxihaai0Wzx3tigqAO8YR05wqzHzk4KjGXZD62z/873KLVNvtdTqfLeoxIYYY/HGRWSTFLMFxf5HTt902lo4Hk0PMzPmJtjktKtWzx2o6N2fkG+m044zEVP0wruiScOXlhM8CzlhC1+C8ALAHqUUjMAXgIQAgCt9ctKqQEANwG0A7CUUr8J4ILW2ocrl/vAlLN99IgXZjnlbC2LlpcJN9SaF9u5c7w492s5aU3BNRb42ho/q67O7u/Z20srL58lmU4znT1TwE00S309/+/4cTtiw+sddYqJfbk3m3R697bLy4xQaWvjjbutrfBNp9B7HOSmY2rV377NG4rp39rdzbHU19vinkzyXLx6lVa9EGhKCrrW+qslXl8AUHtNA/dTznZryw43jMe57alT/N/9Wk2RiO1GWV7mBQzQuj9+nALe1ZV/MTMSyRbvcJhiAXDhtbfXFu+2Nv9N0c14nVzIffiQN80LF1jD/KDWbrEbS6nXUimGHc7N0UDY2uL6S18fZ2AmAurECWkFVyOIy2WvWBat8bt3yytnm07zgpua4sWnFJN+jh3bX7ZgIpEt4NvbfP7QIVr3psVarn81mdxtfRvxb2igaJ86ZQu4CEA2lsUFxakpLm5fueLMuoZSvOHs96Zz7BjHkk7bLfWWl/l7ZyfPz5UVng8mosZvN2ahbETQy2Wv5WzX13nxz87Skmpt5fYjI3tL3LAs3giMG2Vjg883NHAMJ0/yxpB5Q9Ga1namgG9u2q+3tVH8jXgfPixha8VIJOi3XlmhS+3sWe99X4Va6t2+vXtbE+ve3Jwt9KZEgfldqir6DhH0cii3nG0yaYcbhsO8yAYHaUWVG25oxDgzrd6yKCBHjlBMent5IzGikkhwym3Ee33dDkFsbOT/DQ3Z1QHlQi2fcJiZn7GYf/zQmS31TGx7LMaH+d38XFvjz3z+/Pr63cKfewMwbh2v3eBqFBH0YpRTzlZrvj41xbhhy6LoX7rEi78c8YxGs9PqEwk+39bGiJeeHi54NTTw/cNhLsoZ8Y5EuL1Sdhq5sb692MPTLywssCFFQwNrsvgxZttY3aUippLJ3WKf+fvmJs9PYyhkolRh4c+8ATQ2irunwoig56OccraxGP2V09MUVFNlb3SUgl6MZNIOJ1xasgW5uZkLWiapp7nZtqLu3rUF3FhTzc0U7WPH+LOjQzI5neL+fZ4DHR2sYe71iJ6DYko5lKqumE7vtvQzbwAmzNUYJbk0NpZ29TQ3y3m8T0TQM8ktZ3v2LCMEzOKX1na44eIi/+7uZnTL4GDhk9CyeJIbC3x93U7j7+mhFW7iwU3Szscf29NhgDeTjg5ua6xv6efoPJbFErIzM5yNPfmkiEsmpvREqZmfZdm1h/K5euJxWv3xuB1ZlUlDQ2lXj/j5dyGCDtBiNp12gN3lbCMRO9wwFuOJdPIkrfFCJ7aZomaWl1WK0/bTpyngjY22gL//Pl0pxvpuaeHNwoh3e7tMVytNPE5/+doaI5eC3nGpktTV2QJcDK15/RVy9ZiyEvG4XSAu3+eU8vU3NtaEn7+2BT1fOduzZymm6TQjBR49oiArRXfIpUv8mSuusVh2Wn1uedkjR2jdbG3xBJ2ctKel9fUU+hMnbAGXEqbVZWODPT+TSWmCXE3UThGycsJkjZ+/kK9/czM7FyP3c/K5dvL99LHhVJuCrjXdKp99xpOhv58WWXs7reSPPqKYJ5MU93PnKMqZ1kYqRaE3Im7CAhsb6UZpaeG0MRpl2OGjR/bU8vBhfmZm0k4NWA+eZW6OM6TGRuD5591vOCHkZy9+/nwuHvMzGqXb0xhduZgql6VcPh50xdWeoM/Ps+bK1hbF9No1CursLC/qjQ3eoQcH6VLp7qbYmip4mWn1Wtu+7f5+HuBEgq8bKyEUovU9MGALuPj9vIHWLNvw2Wc8Lk89JTOjIFBfT4OqVBZvZj36Qi6frS3+ni+ss6GhvOieKl7vtSPoKyuMWlhbo4A/9RTvxI8eUeTTaVpm4+N0vYRCPJiTk9nlZbXm/x06RDE3C54Ahb+tjYtpmWGDYn17j3SaN/C5OR7vy5d9PdUW9kFmPfpSkWmJRGGrPxajIfj4cWE/f67Q9/XRCHSY4At6OEwhX1zkF3n+PEX59m0udjY08IIeHaVILy+zvvXyMqdmpluP6c5jFnGSSR6cI0f4v0eO0BJ3Ih1cqCyxGP3lpibLyZNuj0jwOsbPX6oefSq129LPvAFEIjQuGxtF0PdEZjnbhgbeEbXmc1ozc/PECYr86ipD1dbX7bZcRsBN3GxdHS14Y3kfOSJlSP3I2hojWdJpxpdX4KISapiGBj5cSugLnqBnlrM1FQ0BWuimDnhzM63vW7dowW9v2+VG29roN29psYW7s1OSdoLAzAxv3M3NwLPPHqz7jyB4kOAIuilne/++3WigpYVi3dBAMd/YoNBHIrTSje/MtHbLtL6DnhlYS5iZ2f37vFlPTEg1SSGQlNPg4lUAXwSwqLUez/O6AvAHAH4RwDaAr2mtbzk90IKYcrY/+QkXNxMJCrkpM2tWslMpukw6Opgqnxs2KAtiwSSV4kzs8WMe9/FxOdaCe0QiDLTo6XHNh/4agG8AeL3A678A4PTO4zqAf77zs7JozezN73+fC5zr63Zd6VCIot7RQct7dDQ78kSss9pge5uLn1tbFPLjx90ekVCLmJIhk5MMaTbJVG4Iutb6LaXUWJFNvgTgda21BvC2UqpTKTWotZ53apBZhMPAn/wJI1fu3uUCZns7v5zBQbpPxsYo4l1dUuu7VllZYQ1zrYHr1/OXOxaEShKP0+icnKRONTczSXF0tGL5Dk740IcATGf8PbPzXGUE/bd+C/jjP6YVPjTEbi3PPMOfFy5Ilp/Ai+ijjxiGev26lBAWqsvKCkV8ft7uF3zpEo3OChuXVV0UVUrdAHADAEZHR/f3Jn/v79mLmtvbdLW8+SatsbY2fmlDQ7bFPjBA14tkAAYfk1/w8CEvomvXJCtXqA6pFKOoJidZBiQUootvbKyqBoUTgj4LYCTj7+Gd53ahtX4FwCsAMDExkadmZhk8/TTw7W9T1D/8kCGKq6uMK45EGJ44PW1XYTt0iP70zk6KvBH6zk7eAOSCDwbJJPDee/RRnjjB2Zq42oRKEw4zKGNmhqLe0cGs46EhV8KcnRD0NwB8XSn1bXAxdKNi/vNMWlsZS3zyJGuHRyJ8rqeHX/LCAq33zU1+0aY924cfMjKmuZlC39VFkTe1Vtrb6XeXjE//EIlw8TMS4cW039mfIJSDZdGdMjlJY7KujgI+NuZ6V6tywha/BeAFAD1KqRkALwEIAYDW+mUA3wVDFu+DYYu/XqnB5qWvD3jhBZbBvXuX1vmJE5xur67SYl9cZBhjImFb5MmkXcp2dtbODm1tpdD39FDku7poybe1UegluchbLC3RMleKN/jubrdHJASVaJTW+NQUFzxbWjgTHBnxTOSc0vm6hVSBiYkJffPmTWffNB5n9Mv0tF23ZXiYvtX1dYYOLS4ywQigO6ajgz/r6ykOCwsU+WiUD8vi66ZLi+nv2d5uC31rq8Q2u8HkJGdnhw/TFSelGASn0Zr5LJOT1A+tOaM3XcZccOsppd7TWk/kfS1Qgm4wLdzW1+lGGR/PngrFYhT2x4/tKop1dRRqk3CUTNJ1s7bG7ZaWbJFPpewSna2tFBRjyWcKfUuL+HErgWXx+D56xON19aq4yARnSSZpGE5O0pXX2MjEtNFR1w2H2hN0gHfSmRla7PE4D8S5c7ujXSyLYUZG4E3DZtOEoq+PYm1ZFPj1dbtt3PIyI22iUbpzTD0YY80fOpQt8EbwpRfo/kkkGNG0sgKcOsVjWus3TdPJJ19N71TKXi86dMgOEmhqku8tH+vrFPHZWV7zXV20xgcHPTMLr01BN6RS9K0/fEgr7swZHqBCBycSsV0zKys8qA0NnF4ZgTc3hXSai64bG7bQr6/zPUzpXYCf1dBgX1ShULbIm4fUjynO5iYXP2MxLn4OD7s9osqhdXYN7kLdd2Kx/M0X6uvtrjqZ56Khrs4WeCPymT9NhdFaIJ1mXfzJSTvjfHiYOuHBvJbaFnTD1hbrnC8u0voeHy+dPZhK0Qo3Ah+L8fnOTrtAfUdHtqVjWfys9XVa9BsbfJjWV/G4bckrRaE3zTKM0Oda9R5ZcHGVxUUuftbXsznJkSNuj2h/WFbx9miZNbTzXZuhUPG2aOZnrgsqlbJdhmZWaX5Go/a5bTA9OHOFPtPS93uAQCRiL3Imk9SFsTG7wY1HEUHP5PFjCnskwiiWixfL94mFw/z/x4/tLkVNTRT3vj7eIPKdCFrz8zIt+XDYbpQRi1HQMy+QUMj+u6kpv9DXit/4wQMmDHV0UMy96LJKp4t3rjc/TWPwXDIbGBdqadbUVDkRtazdIp8p/LHY7htMU1NhC9/MRL1Gvroqg4MUcp9ESImg52JZdMHcvcsDfPIkcPr03i4W0zvUWO/JJE8OE9fe11e63vb2tm3BG6E3F7x5v/p6e+prWdmzgUOHdvvngxRaaVnMG5ie5kV35Ur1983ccEtZ1KnU7v/NbT1WSLD94M/Wmvuaz8I3P3NdP6FQccGvZva2qavy6BHH2txsL3L6zNUpgl6IWIyLpjMzPKgXLjBBYK+YBtJG3MNhPt/SYrtmurvLE6NYLNuKX1+3p8Na8z0aGmyhN1P4zIuptXW3f/7wYX/5RONxLn6urnLd48wZ50QvszlwsQbBhZoDG/90KYu6llxlxudfyMI30WGZ1NcXFnynFm5XV+26KpbFsOOxMV6TfroeMhBBL8XqKsPgNjZoYY+Pl24aW4xoNDssMp3mydvTYwv8XtwG8Xh2hM3GBi8SgxER45dXihfX1pY9TVaKQp/rtvFiE+twmIufiQTw5JMsfVwO+fzT+SzqRKKwf7qYQJufteLqcppksriFn+uOUqq4hW/WnnJJpRilMjnJcykUYvLPsWM0bHyOCHo5aM2p/aef8sQ6dowhcQe1siyLom4E3gixKSRmwiL3KqomTj7TXbO1Zb9uGtoaK0cpnuibm3ZoJsAL4vDhbLdNWxsvFjeEfmGBDSlCISYLdXSUbrxrns+N5ADsxb1Ml0chy9qnFltgSKfzW/iZfvxcMkMyTQjy6qptQJ065VpdlUohgr4Xkkn61j//nJbY2bOcojklbltb2WGRWlO8MsMi93sTSaWyI2s2Nijg5hiHQraf3TTBNqGXm5u8cAz19bv9806HVmaG5cXjbBP32WcU1xMnuE0sxjHmkuufLmRZNzZ6bwYi7A+zcJsp8pEIDbFHj2g4AfY61uHDPP6FXDpeXbgtgQj6ftjcZDTM0hKFbHycd3wnSaWyF1bjcT5/5Eh2WORBsKzdsfLhsO0brq/nZ3R02CUMTFTO5ia3NeMCCsfQmwWuTP90PtdHPv90Os3p8coK9/v8eV5w+RYPjXD78EIUHCQWo4ibZvAtLVzg7O+3QzTzuXZy10RM2HA+l47x43sMEfSDsLBAYd/eZqTFhQuVSf3VmmJrXDPr63y+uTk7LNIJ/62Jlc902YTD9qJVXR1FuqODMffG37+9bVvzJuzSYCzhQv7pxsb8/miAbq7tbeCJJzgjEoR8OFFXxeSDFIrHL5SAVcjCb26u+gxQBP2gWBZjoe/d40l06hQflfTLxeN2pcjFRbvejJlO9vc7Wzg/s2FIpsvGnOBK2a4XI/RNTdkiDxQOy8vnn15fB959l59x9SrzAgQhl3x1VUZHuc7ltHGVu3CbK/z5Fm4z/fi5Fn6hhdsDIILuFLEYE1xmZ3mgLlwoPwLjIFhWdlikEc/W1uywyEos6kWj2e4ak/VqaG2luBu3TUdHee6Q2Vng/fcp9k8/7ckUa8FlNjbsuirpNF2RY2O85txawE6nC1v429v5M3xNo51MkT9yZN/uVBF0p1ldZc/KcJhCOj5eXUHa3s4Oi7QszhYyF1YrmSwRi2XHyW9sZC+omrLEmUKf6WP/7DPOdrq6mPlZS/HaQnEsy66rsrbG89o0jzjoelI1sKzCCViZJblPn2YU3T4QQa8EWjPz7M4dTtOOHaP/t9rilE5nh0UaYW1vt8X9yJHK+/kSiWxXzcZGdnhkczOt+bk5nuBnz1LMJVRQAHhOTE7StZJI+Kauyp4xQQMmnHYfHFjQlVIvAvgDAPUAvqm1/ic5rx8D8CqAXgCrAP4TrfVMsff0vaAbkklanJOTXLA8d47i7lao3Oam7ZpZXeUJ1NiYbb1X6wJJpWxxf/wYeOcdjmlkhP7yUGi3u0ZqyNcOWvM8nZzkT6V4XoyNOR9RFiAOJOhKqXoAdwH8PIAZAO8C+KrW+nbGNn8M4Dta6z9USv37AH5da/1rxd43MIJuCIcZDbO8TOt4fNz9Yj/JZHZYZCLBiyYzLLIarqLVVS5+as3Mz6am3bHyJpysoSFb4Ds6aK2JyAeHeNxe5DR1Vcwip8/qqrjBQQX9WQC/o7X+D3b+/q8BQGv9exnbfALgRa31tFJKgY2iiypF4ATdMD9PYY9GuXhz4YI3qgOaNnzGNWPa8DU325Z7b6/zkTvT0yywdegQFz/zpV5nlhw2Ih8O2wlFdXW7Rb6tTdw1fiOAdVX2TCLB6+/w4X2XgC4m6OUENQ8BmM74ewbA9ZxtPgDwH4JumV8G0KaU6tZar+QM5AaAGwAwGtTO7IODFEcT5vj4MUMcT550N/3YWOZHjtB/bdrwLS4yiuDRo+w2fH19BwuL1Jrx5Q8e8MKdmCjs6qmr40whc7agNUU+05KfmaEgmP8xsfLm0d4eqBTvQJBbV6WhgZb42Fgg6qqURTTKfJb5edsNeuJERWr6l2Ohfxm0vv/uzt+/BuC61vrrGdscBfANAMcBvAXgVwCMa63XC71vYC30TKJRhjnOzdkdwgcH3R7VbiyLJ5pxzZiaMIcP266Zrq7yrahkkvVYFhd54V686IwFZmLlcxdfTWywUhxzrsgHaVHNL2xt2YucqRSPw/HjnLXWQnGzrS0K+Py8PRtua+MaweDggSJ2DmqhzwIYyfh7eOe5n6K1ngMtdCilDgP4lWJiXjMcOgRcu0ZR+/hjloPt6aF/vVSt9GpSV8dx9fRQfCMR2zUzOWm37+vttQW+0Ap9JMJKiZEIcOkS990pTMXI1tbs+H8TK28ey8u05g2trbtdNhIq6TyWZZ8zy8s8r0zziK4ut0dXWUym9/w8rXFjFB05wlIWg4POJgIWoBwLvQFcFP05UMjfBfCrWutPMrbpAbCqtbaUUv8DgLTW+reLvW9NWOiZaE23xp07tFjGxuj68Lr1aNrwGYE3Fe86OuyMVdOGb3mZNy2ALhY3IxXi8d2WfGbJYRMrn/mQBbn9YeqqTE3x90OHeH6PjHiyFopjmJmtEfFYjNdBTw8t8YGBipxTB7LQtdYppdTXAfwlGLb4qtb6E6XU7wK4qbV+A8ALAH5PKaVBl8t/7tjog4JSdpabCXOcnWWY4+iod6M4GhrskxOw2/AtLnKN4O5dWrtdXXy+tZWLn1WwRoqS2RrQkEzuFvmFhez/yRX5StTtCQqmrsrCAg2Wvj7W4+nr8+75fFDSaUaOzc/zfE8m7aS+wUEaOC4aaZJY5BbhMLNNV1cpHOPj/puWZrbhW17m9PLJJ70/68iknJLDuSLvxaYg1SKZtBent7YqW1fFKySTPMfn53m+p9M8L/r7KeKViA4rgmSKepm5OYY5xmJMcb5wQab+bmNZu0W+UMnhzFj5IIfehcPsEeCluiqVJBbjzGNhgcaK1rwuzWy1UrWTyuCgi6JCJTl6lHf6e/cY4rewwDoPJ08G80LxA3V1zGDt7LSfM7HymSI/NZUdK28qUWZG2Pj5GPq9rspeiURsf/jaGp9rbeW1ODDA88HjMzOx0L3E9jat9YUFTl8vXpSSsl7GNALJ9cvnlhzOFXmvh+1tb9uLnIkERc0scvrJnVYOmZEppoppRwddKQMD3opG20EsdL/Q0sKCVUtLFPZ336V/bny8dpIw/ISJez98mJarITdWfnGR8diG3Fj5cksOVxKted6Z5hFKceZ4/Hiw6qpozXUrk+gTjXJfTXP4gQFvZHbvE7HQvYpl8eK6e5cLd8ePA2fOuH/hC/vDlBzOfGSWHG5p2S3y1Qj5SyRoiT96xBtRUxMXOINUV8WyeLMyPvFEgq6wzMgUH+UliIXuR+rqmB48PMwU+ocP7TDHkRHP+/KEHEwnp/5++7l8JYfn57P/J1fknbIe19ZoMMzNUfC6u5kAMzDgb7+/IZXiTGNhwe741dDA739ggKGVXnd97YPg7VHQaGwELl+2s00/+IDW1Ph4RWpBCFXElDXu7bWfSyZ3R9gsLtphlI2N+WPly7nBp9N2XZWNDbuuyrFjnvQV75l43A4vNI1fmproDhsYoOsoCDerIoig+4WODuC553hB3r4N/OAHtN7Pnw/O1FigS627O7v0cjq9W+QfPiy/5PDWFo2A6WneMNrbmQA0NOR/K3V7O7vwFcAb3PHjFPFqNHfxED4/mjXI0BCnjffvZ4c5njgReOujZqmvtytlGiyLURmZIv/okR1GWV9vh02urASrrko4bIt4OMzn2ttZSmNgoKb704qg+xHTGWlkhNb6p59yYevixWwfrRBcMmvEG/KVHE4k7PISfq2rojV9/kbETU2eri47tDeoWap7RATdz7S2MsxxcZFhjj/+MRd7xsfdr6UiVB8T997WRnecn7Es+sFNZEo8blcFPXWKIu7XG1QFEUEPAn19PNEnJ1n468037TBHv/tIhdohlaJxsrDAxU0TmdLXZ0emSNhuUeRqDwomzHFoiCV6HzxgEaXz52mt1dDCkOAjEgnbCl9aomXe2MiSGAMDjACStaGyEUEPGk1NDHM8doxhju+/b4c5ZtYmEQS3yNeSzdRQHxigb1wMkH0hgh5UOjuzwxy//30uop4/L75HofpsbtointmS7dSpA7dkE2xE0IOMUnS3DAywhMDnn/OCOnOGPnaZygqVQmtgfd0W8UiEz1e5JVutUZagK6VeBPAHYMeib2qt/0nO66MA/hBA5842/0hr/V1nhyrsm4YG1lk3bpjbt+0wx8yOPoJwECyLMe/GJ57Zku3EiYq1ZBNsSgq6UqoewP8G4OcBzAB4Vyn1htb6dsZm/x2AP9Ja/3Ol1AUA3wUwVoHxCgehtRW4fp0RBJ98ArzzDi+yCxfEWhL2h8dbstUa5VjoTwO4r7V+CABKqW8D+BKATEHXAEx6VgeAOScHKThMfz8vuIcP2VjjzTdZxP/UKQlzFErjsZZsgk05V+8QgIxizpgBcD1nm98B8G+VUn8fQCuAv57vjZRSNwDcAIDR0dG9jlVwkro6Crip5njvHmt9XLiQXdtbEIDCLdlGRijiXV2yJuMBnDLHvgrgNa31/6iUehbAv1BKjWutrcyNtNavAHgFYD10hz5bOAjNzcCVK3Y1x1u3mKA0Pi6RB7VOvpZshw9zNmciUyS80FOUI+izAEYy/h7eeS6TvwPgRQDQWv9IKdUMoAfAohODFKrAkSPA88/TSr9zB3jrLdb/OH/eV8X/hQOSryVbZyfrwQwOSucsj1OOoL8L4LRS6jgo5F8B8Ks520wB+DkArymlzgNoBrDk5ECFKqAURXxwcHeY49iYTKmDiGnJZkTctGTr7mZUlM9bstUaJQVda51SSn0dwF+CIYmvaq0/UUr9LoCbWus3APyXAP4PpdRvgQukX9Nu9bYTDk4oxJBGE+b4ySd2mGNmMwbBn5iWbCYyJbMl29mzvmvJJthIT1GhNAsLFPXtbVpsFy9KuVK/kUxmF75Kp+2WbCYyRSKcfIH0FBUOhql09+ABo2G+9z1GyJw6JeFp1cay6CaxrN0PrSnUma+buimZLdmGhyni3d3iRgsYIuhCedTVsTOSaapx967thjl61O3ROUOmMOYTzEKimW+bQq8d9D32Q2trzbZkqzVE0IW90dwMXL1qhzm+954d5lis9ZdXBLHY65Wgrm73Q6ndzzU0lN5mL6+bbRobJQu4hhBBF/ZHVxfwMz/DMMdPP2WYY0tLYdGsBG6LZX194W3N84JQRUTQhf2TGeb44AEXTcsRy2JCKGIpCPtGBF04OKEQE08EQXAVMXMEQRACggi6IAhCQBBBFwRBCAgi6IIgCAFBBF0QBCEgiKALgiAEBBF0QRCEgCCCLgiCEBBcK5+rlFoC8Gif/94DYNnB4biJ7Is3Ccq+BGU/ANkXwzGtdd7GBK4J+kFQSt0sVA/Yb8i+eJOg7EtQ9gOQfSkHcbkIgiAEBBF0QRCEgOBXQX/F7QE4iOyLNwnKvgRlPwDZl5L40ocuCIIg7MavFrogCIKQgwi6IAhCQPC0oCulXlRKfaaUuq+U+kd5Xm9SSv3LndffUUqNuTDMsihjX76mlFpSSr2/8/i7boyzFEqpV5VSi0qpjwu8rpRS/8vOfn6olLpa7TGWSxn78oJSaiPjmPx2tcdYDkqpEaXU95RSt5VSnyil/os82/jiuJS5L345Ls1KqR8rpT7Y2Zd/nGcbZzVMa+3JB4B6AA8AnADQCOADABdytvnPALy88/tXAPxLt8d9gH35GoBvuD3WMvblZwFcBfBxgdd/EcBfAFAAngHwjttjPsC+vADgO26Ps4z9GARwdef3NgB385xfvjguZe6LX46LAnB45/cQgHcAPJOzjaMa5mUL/WkA97XWD7XWCQDfBvClnG2+BOAPd37/VwB+TimlqjjGcilnX3yB1votAKtFNvkSgNc1eRtAp1JqsDqj2xtl7Isv0FrPa61v7fy+CeBTAEM5m/niuJS5L75g57ve2vkztPPIjUJxVMO8LOhDAKYz/p7B7gP702201ikAGwC6qzK6vVHOvgDAr+xMh/+VUmqkOkNznHL31S88uzNl/gul1EW3B1OKnSn7FdAazMR3x6XIvgA+OS5KqXql1PsAFgH8O611wePihIZ5WdBrjT8HMKa1fgLAv4N91xbc4xZYN+MygP8VwJ+5O5ziKKUOA/gTAL+ptQ67PZ6DUGJffHNctNZprfWTAIYBPK2UGq/k53lZ0GcBZFqpwzvP5d1GKdUAoAPASlVGtzdK7ovWekVrHd/585sArlVpbE5TznHzBVrrsJkya62/CyCklOpxeVh5UUqFQAH8v7XWf5pnE98cl1L74qfjYtBarwP4HoAXc15yVMO8LOjvAjitlDqulGoEFwzeyNnmDQD/6c7vXwbw/+md1QWPUXJfcvyZvwT6Dv3IGwD+9k5UxTMANrTW824Paj8opQaMP1Mp9TR4vXjOYNgZ4/8J4FOt9f9UYDNfHJdy9sVHx6VXKdW58/shAD8P4E7OZo5qWMN+/7HSaK1TSqmvA/hLMErkVa31J0qp3wVwU2v9Bnjg/4VS6j64uPUV90ZcmDL35R8opX4JQArcl6+5NuAiKKW+BUYZ9CilZgC8BC72QGv9MoDvghEV9wFsA/h1d0ZamjL25csAfkMplQIQBfAVjxoMzwH4NQAf7fhrAeC/ATAK+O64lLMvfjkugwD+UClVD950/khr/Z1Kapik/guCIAQEL7tcBEEQhD0ggi4IghAQRNAFQRACggi6IAhCQBBBFwRBCAgi6IIgCAFBBF0QBCEg/P/Spz8kabH7swAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.09</td>\n",
       "      <td>1.08</td>\n",
       "      <td>1.34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.16</td>\n",
       "      <td>1.26</td>\n",
       "      <td>1.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.22</td>\n",
       "      <td>1.07</td>\n",
       "      <td>1.03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.93</td>\n",
       "      <td>0.97</td>\n",
       "      <td>0.92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.11</td>\n",
       "      <td>1.56</td>\n",
       "      <td>1.52</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.76</td>\n",
       "      <td>0.77</td>\n",
       "      <td>0.90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.92</td>\n",
       "      <td>0.84</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.88</td>\n",
       "      <td>1.22</td>\n",
       "      <td>1.34</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     0     1     2     3\n",
       "1  1.0  1.09  1.08  1.34\n",
       "2  1.0  1.16  1.26  1.54\n",
       "3  1.0  1.22  1.07  1.03\n",
       "4  1.0  0.93  0.97  0.92\n",
       "5  1.0  1.11  1.56  1.52\n",
       "6  1.0  0.76  0.77  0.90\n",
       "7  1.0  0.92  0.84  1.01\n",
       "8  1.0  0.88  1.22  1.34"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "paths = [\n",
    "    {0:1.00, 1:1.09, 2:1.08, 3:1.34},\n",
    "    {0:1.00, 1:1.16, 2:1.26, 3:1.54},\n",
    "    {0:1.00, 1:1.22, 2:1.07, 3:1.03},\n",
    "    {0:1.00, 1:0.93, 2:0.97, 3:0.92},\n",
    "    {0:1.00, 1:1.11, 2:1.56, 3:1.52},\n",
    "    {0:1.00, 1:0.76, 2:0.77, 3:0.90},\n",
    "    {0:1.00, 1:0.92, 2:0.84, 3:1.01},\n",
    "    {0:1.00, 1:0.88, 2:1.22, 3:1.34},\n",
    "]\n",
    "\n",
    "df = pd.DataFrame(paths, index=range(1, len(paths)+1))\n",
    "\n",
    "mu = 0.06\n",
    "k = 1.10\n",
    "\n",
    "df.transpose().plot(color=\"red\", alpha=0.3)\n",
    "plt.legend([])\n",
    "plt.plot([0,N], [k, k], label=\"strike price\")\n",
    "plt.show()\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generated stock price evolution and strike price. Now for each of these $M$ paths we compute the option value\n",
    "\n",
    "### last timestamp\n",
    "\n",
    "we can easily compute the option value\n",
    "\n",
    "We compute the discounted payoff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> If the put is in the money at time 2, the optionholder must then decide whether to exercise the option immediately or continue the option’s life until the final expiration date at time 3\n",
    "\n",
    "We need to compare what is worth more: exercising the option at expiry, or exercising it prematurely. Since we talk about exercise value (not extrinsic value), we only consider ITM options, since exercising an OTM option does not make sense (but selling does, however that is not what we are doing).\n",
    "\n",
    "By ITM options we mean options that are ITM at time 2, since that is where we face the decision of early exercise or continuation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    0.000000\n",
       "2    0.000000\n",
       "3    0.065924\n",
       "4    0.169518\n",
       "5    0.000000\n",
       "6    0.188353\n",
       "7    0.084759\n",
       "8    0.000000\n",
       "Name: 3, dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = (k - df[N]).map(lambda v: max(v, 0))  # put payoff \n",
    "\n",
    "# discount the payoff\n",
    "discount = np.exp(-mu * 1)\n",
    "Y = Y * discount\n",
    "\n",
    "Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "X is the stock value at current time step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    1.08\n",
       "3    1.07\n",
       "4    0.97\n",
       "6    0.77\n",
       "7    0.84\n",
       "Name: 2, dtype: float64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = df[N-1]\n",
    "\n",
    "ITM = X < k\n",
    "\n",
    "X = X[ITM]\n",
    "\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    0.000000\n",
       "3    0.065924\n",
       "4    0.169518\n",
       "6    0.188353\n",
       "7    0.084759\n",
       "Name: 3, dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = Y[ITM]\n",
    "\n",
    "Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have X and Y\n",
    "\n",
    "X is the stock price at time 2\n",
    "\n",
    "Y is the option price at time 3, discounted to time 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fd7b2715bb0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVo0lEQVR4nO3db4xd9X3n8fenAzh+sMQEphLYJHYWl9bdSNBeaLRRqDZRYlO12JuSxDSbkC4qm+2yT6Kg2Iq6XXlbNREPWEVi27jNXxRCWESIpWzkzS5JVooK9TUmGCdyMzhp8ICUCWDaVbyA4bsP7jG53Iw9Z/54/vi8X9KRz/n9OfP7cdD9zD3n3vmlqpAkdc8vLfUAJElLwwCQpI4yACSpowwASeooA0CSOuqcpR7AbFx00UW1fv36pR6GJK0o+/fv/2lVjY+Wr6gAWL9+Pf1+f6mHIUkrSpJ/mK7cW0CS1FEGgCR1lAEgSR1lAEhSRxkAktRRK+pTQHNx/4FJbtt7mCePHeeSNau5dfPlbLty7VIPS5KW3FkdAPcfmGTnfQc5/uJLAEweO87O+w4CGAKSOu+svgV0297Dr7z4n3T8xZe4be/hJRqRJC0fZ3UAPHns+KzKJalLzuoAuGTN6lmVS1KXnNUBcOvmy1l97tirylafO8atmy9fohFJ0vJxVj8EPvmg108BSdIvOqsDAAYh4Au+JP2is/oWkCTp1AwASeooA0CSOsoAkKSOMgAkqaMMAEnqqFYBkGRLksNJJpLsmKb+miQPJzmR5Pqh8n+V5JGh7f8l2dbUfS7JD4fqrlioSUmSZjbj9wCSjAF3AO8AjgL7kuypqu8NNfsx8EHgI8N9q+qbwBXNeV4HTAD/c6jJrVV17zzGL0maozZfBLsamKiqIwBJ7ga2Aq8EQFX9qKl7+TTnuR74elX9bM6jlSQtmDa3gNYCTwwdH23KZms78KWRsj9P8miS25Osmq5TkpuT9JP0p6am5vBjJUnTWZSHwEkuBt4E7B0q3gn8KnAV8Drgo9P1rardVdWrqt74+PgZH6skdUWbAJgELh06XteUzcZ7gK9U1YsnC6rqqRp4Hvgsg1tNkqRF0uYZwD5gY5INDF74twN/MMufcwOD3/hfkeTiqnoqSYBtwGOzPKekWXKNbA2b8R1AVZ0AbmFw++b7wD1VdSjJriTXASS5KslR4N3Ap5IcOtk/yXoG7yC+PXLqLyY5CBwELgL+bAHmI+kUTq6RPXnsOMXP18i+/8Bs39DrbJGqWuoxtNbr9arf7y/1MKQV6S0ff4DJaZZDXbtmNd/Z8bYlGJEWS5L9VdUbLfebwFJHuEa2RhkAUke4RrZGGQBSR7hGtkad9UtCShpwjWyNMgCkDnGNbA3zFpAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRrQIgyZYkh5NMJNkxTf01SR5OciLJ9SN1LyV5pNn2DJVvSPJQc84vJzlv/tORJLU1YwAkGQPuAK4FNgE3JNk00uzHwAeBu6Y5xfGquqLZrhsq/wRwe1VdBjwL3DSH8UuS5qjNO4CrgYmqOlJVLwB3A1uHG1TVj6rqUeDlNj80SYC3Afc2RZ8HtrUdtCRp/toEwFrgiaHjo01ZW69J0k/yYJJtTdmFwLGqOjHTOZPc3PTvT01NzeLHSpJOZzEWhHlDVU0meSPwQJKDwHNtO1fVbmA3QK/XqzM0RknqnDbvACaBS4eO1zVlrVTVZPPvEeBbwJXA08CaJCcDaFbnlCTNX5sA2AdsbD61cx6wHdgzQx8AklyQZFWzfxHwFuB7VVXAN4GTnxi6EfjqbAcvSZq7GQOguU9/C7AX+D5wT1UdSrIryXUASa5KchR4N/CpJIea7r8G9JN8l8EL/ser6ntN3UeBDyeZYPBM4NMLOTFJ0ull8Mv4ytDr9arf7y/1MCRpRUmyv6p6o+V+E1iSOsoAkKSOMgAkqaMMAEnqKANAkjrKAJCkjjIAJKmjDABJ6igDQJI6ygCQpI4yACSpowwASeooA0CSOsoAkKSOMgAkqaMMAEnqKANAkjrKAJCkjmoVAEm2JDmcZCLJjmnqr0nycJITSa4fKr8iyd8mOZTk0STvHar7XJIfJnmk2a5YkBlJklo5Z6YGScaAO4B3AEeBfUn2DC3uDvBj4IPAR0a6/wz4QFX9IMklwP4ke6vqWFN/a1XdO885SJLmYMYAAK4GJqrqCECSu4GtwCsBUFU/aupeHu5YVX8/tP9kkp8A48Cx+Q5ckjQ/bW4BrQWeGDo+2pTNSpKrgfOAx4eK/7y5NXR7klWn6Hdzkn6S/tTU1Gx/rCTpFBblIXCSi4E7gT+sqpPvEnYCvwpcBbwO+Oh0fatqd1X1qqo3Pj6+GMOVpE5oEwCTwKVDx+uaslaSnA98DfhYVT14sryqnqqB54HPMrjVJElaJG0CYB+wMcmGJOcB24E9bU7etP8K8IXRh73NuwKSBNgGPDaLcUuS5mnGAKiqE8AtwF7g+8A9VXUoya4k1wEkuSrJUeDdwKeSHGq6vwe4BvjgNB/3/GKSg8BB4CLgzxZyYpKk00tVLfUYWuv1etXv95d6GJK0oiTZX1W90XK/CSxJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR3VKgCSbElyOMlEkh3T1F+T5OEkJ5JcP1J3Y5IfNNuNQ+W/meRgc85PNmsDS5IWyYwBkGQMuAO4FtgE3JBk00izHwMfBO4a6fs64E+B3wKuBv40yQVN9V8CfwRsbLYtc56FJGnW2rwDuBqYqKojVfUCcDewdbhBVf2oqh4FXh7puxn4RlU9U1XPAt8AtiS5GDi/qh6swaLEXwC2zXMukqRZaBMAa4Enho6PNmVtnKrv2mZ/xnMmuTlJP0l/amqq5Y+VJM1k2T8ErqrdVdWrqt74+PhSD0eSzhptAmASuHToeF1T1sap+k42+3M5pyRpAbQJgH3AxiQbkpwHbAf2tDz/XuCdSS5oHv6+E9hbVU8B/5jkzc2nfz4AfHUO45ckzdGMAVBVJ4BbGLyYfx+4p6oOJdmV5DqAJFclOQq8G/hUkkNN32eA/8IgRPYBu5oygD8G/gaYAB4Hvr6gM5MknVYGH8JZGXq9XvX7/aUehiStKEn2V1VvtHzZPwSWJJ0ZBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHXXOUg9Ai+/+A5PctvcwTx47ziVrVnPr5svZdmXbv+8n6WxhAHTM/Qcm2XnfQY6/+BIAk8eOs/O+gwCGgNQx3gLqmNv2Hn7lxf+k4y++xG17Dy/RiCQtFQOgY548dnxW5ZLOXgZAx1yyZvWsyiWdvQyAjrl18+WsPnfsVWWrzx3j1s2XL9GIJC0VHwJ3zMkHvX4KSJIB0EHbrlzrC74kbwFJUlcZAJLUUa0CIMmWJIeTTCTZMU39qiRfbuofSrK+KX9fkkeGtpeTXNHUfas558m6X17IiUmSTm/GAEgyBtwBXAtsAm5Ismmk2U3As1V1GXA78AmAqvpiVV1RVVcA7wd+WFWPDPV738n6qvrJvGcjSWqtzTuAq4GJqjpSVS8AdwNbR9psBT7f7N8LvD1JRtrc0PSVJC0DbQJgLfDE0PHRpmzaNlV1AngOuHCkzXuBL42Ufba5/fMn0wQGAEluTtJP0p+ammoxXElSG4vyEDjJbwE/q6rHhorfV1VvAt7abO+frm9V7a6qXlX1xsfHF2G0ktQNbQJgErh06HhdUzZtmyTnAK8Fnh6q387Ib/9VNdn8+0/AXQxuNUmSFkmbANgHbEyyIcl5DF7M94y02QPc2OxfDzxQVQWQ5JeA9zB0/z/JOUkuavbPBX4XeAxJ0qKZ8ZvAVXUiyS3AXmAM+ExVHUqyC+hX1R7g08CdSSaAZxiExEnXAE9U1ZGhslXA3ubFfwz4X8BfL8iMJEmtpPlFfUXo9XrV7/eXehiStKIk2V9VvdFyvwksSR3lH4OTpCWwHNbmNgAkaZEtl7W5vQUkSYtsuazNbQBI0iJbLmtzGwCStMiWy9rcBoAkLbLlsja3D4ElaZEtl7W5DQBJWgLLYW1ubwFJUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR3VKgCSbElyOMlEkh3T1K9K8uWm/qEk65vy9UmOJ3mk2f5qqM9vJjnY9PlkkizYrCRJM5oxAJKMAXcA1wKbgBuSbBppdhPwbFVdBtwOfGKo7vGquqLZPjRU/pfAHwEbm23L3KchSZqtNu8ArgYmqupIVb0A3A1sHWmzFfh8s38v8PbT/Uaf5GLg/Kp6sAaLEn8B2DbbwUuS5q5NAKwFnhg6PtqUTdumqk4AzwEXNnUbkhxI8u0kbx1qf3SGcwKQ5OYk/ST9qampFsOVJLVxph8CPwW8vqquBD4M3JXk/NmcoKp2V1Wvqnrj4+NnZJCS1EVtAmASuHToeF1TNm2bJOcArwWerqrnq+ppgKraDzwO/ErTft0M55QknUFtAmAfsDHJhiTnAduBPSNt9gA3NvvXAw9UVSUZbx4ik+SNDB72Hqmqp4B/TPLm5lnBB4CvLsB8JEktzbgeQFWdSHILsBcYAz5TVYeS7AL6VbUH+DRwZ5IJ4BkGIQFwDbAryYvAy8CHquqZpu6Pgc8Bq4GvN5skaZFk8CGclaHX61W/31/qYUjSipJkf1X1Rsv9JrAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHdUqAJJsSXI4yUSSHdPUr0ry5ab+oSTrm/J3JNmf5GDz79uG+nyrOecjzfbLCzYrSdKMZlwUPskYcAfwDuAosC/Jnqr63lCzm4Bnq+qyJNuBTwDvBX4K/F5VPZnkXzBYWH7tUL/3VZWL/ErSEmjzDuBqYKKqjlTVC8DdwNaRNluBzzf79wJvT5KqOlBVTzblh4DVSVYtxMAlSfPTJgDWAk8MHR/l1b/Fv6pNVZ0AngMuHGnz+8DDVfX8UNlnm9s/f5Ik0/3wJDcn6SfpT01NtRiuJKmNRXkInOTXGdwW+ndDxe+rqjcBb22290/Xt6p2V1Wvqnrj4+NnfrCS1BFtAmASuHToeF1TNm2bJOcArwWebo7XAV8BPlBVj5/sUFWTzb//BNzF4FaTJGmRtAmAfcDGJBuSnAdsB/aMtNkD3NjsXw88UFWVZA3wNWBHVX3nZOMk5yS5qNk/F/hd4LF5zUSSNCszBkBzT/8WBp/g+T5wT1UdSrIryXVNs08DFyaZAD4MnPyo6C3AZcB/Gvm45ypgb5JHgUcYvIP46wWclyRpBqmqpR5Da71er/p9PzUqSbORZH9V9UbL/SawJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1VKsASLIlyeEkE0l2TFO/KsmXm/qHkqwfqtvZlB9OsrntOSWp6+4/MMlbPv4AG3Z8jbd8/AHuPzC5oOefMQCSjAF3ANcCm4AbkmwaaXYT8GxVXQbcDnyi6buJwSLyvw5sAf5bkrGW55Skzrr/wCQ77zvI5LHjFDB57Dg77zu4oCHQ5h3A1cBEVR2pqheAu4GtI222Ap9v9u8F3p4kTfndVfV8Vf0QmGjO1+acktRZt+09zPEXX3pV2fEXX+K2vYcX7Ge0CYC1wBNDx0ebsmnbVNUJ4DngwtP0bXNOAJLcnKSfpD81NdViuJK08j157Pisyudi2T8ErqrdVdWrqt74+PhSD0eSFsUla1bPqnwu2gTAJHDp0PG6pmzaNknOAV4LPH2avm3OKUmddevmy1l97tirylafO8atmy9fsJ/RJgD2ARuTbEhyHoOHuntG2uwBbmz2rwceqKpqyrc3nxLaAGwE/q7lOSWps7ZduZa/eNebWLtmNQHWrlnNX7zrTWy7ctq75XNyzkwNqupEkluAvcAY8JmqOpRkF9Cvqj3Ap4E7k0wAzzB4Qadpdw/wPeAE8B+q6iWA6c65YLOSpLPAtivXLugL/qgMflFfGXq9XvX7/aUehiStKEn2V1VvtHzZPwSWJJ0ZBoAkdZQBIEkdZQBIUketqIfASaaAf1jqcczCRcBPl3oQC+BsmIdzWB6cw9J4Q1X9wjdpV1QArDRJ+tM9eV9pzoZ5OIflwTksL94CkqSOMgAkqaMMgDNr91IPYIGcDfNwDsuDc1hGfAYgSR3lOwBJ6igDQJI6ygCYo5kWtU/y+iTfTHIgyaNJfmeobmfT73CSzYs78leNcU5zSLI+yfEkjzTbXy3+6F8Z40xzeEOS/92M/1tJ1g3V3ZjkB81242jfxTLPObw0dB2W7E+qJ/lMkp8keewU9UnyyWaOjyb5jaG65XId5jOHZXEdZq2q3Ga5MfgT1o8DbwTOA74LbBppsxv4983+JuBHQ/vfBVYBG5rzjK2wOawHHlsh1+G/Azc2+28D7mz2Xwccaf69oNm/YCXNoTn+v0t9HZpxXAP8xqn+vwB+B/g6EODNwEPL6TrMZw7L6TrMdvMdwNy0WdS+gPOb/dcCTzb7W4G7q+r5qvohMNGcb7HNZw7LRZs5bAIeaPa/OVS/GfhGVT1TVc8C3wC2LMKYR81nDstGVf0fBmuBnMpW4As18CCwJsnFLJ/rMJ85rFgGwNy0WdT+PwP/JslR4H8A/3EWfRfDfOYAsKG5NfTtJG89oyM9tTZz+C7wrmb/XwP/LMmFLfsuhvnMAeA1SfpJHkyy7YyOdH5ONc/lch3aON1YV8p1eBUD4My5AfhcVa1j8NbxziQr7b/3qebwFPD6qroS+DBwV5LzT3OepfQR4LeTHAB+m8Ha0y8t7ZBm7XRzeEMN/izBHwD/Nck/X6Ixdt2KvA4r7QVpuWizqP1NwD0AVfW3wGsY/BGpNn0Xw5zn0Ny+erop38/gHvavnPER/6IZ51BVT1bVu5qw+lhTdqxN30UynzlQVZPNv0eAbwFXnvkhz8mp5rlcrkMbpxzrCroOr2IAzE2bRe1/DLwdIMmvMXjxnGrabU+yKskGYCPwd4s28p+b8xySjCcZa8rfyGAORxZt5D834xySXDT0zmsn8Jlmfy/wziQXJLkAeGdTttjmPIdm7KtOtgHewmD97eVoD/CB5pM0bwaeq6qnWD7XoY1p57DCrsOrLfVT6JW6Mbgl8vcMfvv9WFO2C7iu2d8EfIfB/dtHgHcO9f1Y0+8wcO1KmwPw+8Chpuxh4PeW8RyuB37QtPkbYNVQ33/L4CH8BPCHK20OwL8EDjbX5yBw0xLO4UsMbg2+yODe+E3Ah4APNfUB7mjmeBDoLcPrMKc5LKfrMNvNPwUhSR3lLSBJ6igDQJI6ygCQpI4yACSpowwASeooA0CSOsoAkKSO+v+eLEc8sVV9bQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we plot X against Y, we see the values of the option as a function of the stock price. We now want to find a solution that we can use the find the value of the option for any stock price. We find such a function by approximating it with basis functions\n",
    "\n",
    "The most simple basis functions are probably $1, x, x^2, x^3, ...$\n",
    "\n",
    "Others are like $1, cos(x), sin(x), cos(2x), sin(2x), ...$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.471328</td>\n",
       "      <td>0.881958</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0.480124</td>\n",
       "      <td>0.877201</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>0.565300</td>\n",
       "      <td>0.824886</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1</td>\n",
       "      <td>0.717911</td>\n",
       "      <td>0.696135</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1</td>\n",
       "      <td>0.667463</td>\n",
       "      <td>0.744643</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   0         1         2\n",
       "1  1  0.471328  0.881958\n",
       "3  1  0.480124  0.877201\n",
       "4  1  0.565300  0.824886\n",
       "6  1  0.717911  0.696135\n",
       "7  1  0.667463  0.744643"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poly = pd.DataFrame(index=X.index)\n",
    "\n",
    "poly[0] = 1\n",
    "poly[1] = np.cos(X)\n",
    "poly[2] = np.sin(X)\n",
    "# poly[3] = np.cos(2*X)\n",
    "# poly[4] = np.sin(2*X)\n",
    "\n",
    "poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sm.OLS(Y, poly)\n",
    "res = model.fit()\n",
    "coef = res.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "continuation = (poly * coef).sum(axis=1)\n",
    "\n",
    "exercise = (k - df[N-1][ITM])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd7b28361c0>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAI/CAYAAABAoBw9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABGa0lEQVR4nO3dd3zUReLG8WdSIAJKUUQEKYogLQEJqChFqgKi2BWPoqee9cRyqChylrOgB4oellOxcD8R7lRUzoIU4UQlNEEQRQWkKEgJBAikzO+PSaGE1N2dLZ/365VX9ru7yT7ZxfhkZna+xlorAAAABEec7wAAAADRjLIFAAAQRJQtAACAIKJsAQAABBFlCwAAIIgoWwAAAEGU4DvA4RxzzDG2UaNGvmMAAACUaMGCBb9ba2sXdVvYlq1GjRopLS3NdwwAAIASGWPWHO42phEBAACCiLIFAAAQRJQtAACAIArbNVsAgOiTlZWldevWKTMz03cUoFySkpJUv359JSYmlvprKFsAgJBZt26djjzySDVq1EjGGN9xgDKx1mrLli1at26dGjduXOqvYxoRABAymZmZOvrooylaiEjGGB199NFlHpmlbAEAQoqihUhWnn+/lC0AAIAgomwBAGKGtVZnnXWW/vvf/xZcN3nyZJ1zzjmH3K9bt27asWPHId9j1KhRevLJJyVJI0eO1PTp04MbuhiLFy/WtGnTyvx1Xbt29b5xeCCeu1mzZqlfv37F3qe8z9HmzZsP+XdRXpQtAEDMMMbo+eef1+23367MzExlZGTo3nvv1XPPPXfA/aZNm6aUlBQdddRRxX6/Bx98UD169Ahm5GKVt0iEirVWubm5Rd4WqueuPM9Rdna2ateurbp16+p///tfhTNQtgAA4WviRKlRIykuzn2eOLHC37JVq1Y677zz9Pjjj+vBBx/UoEGDdNJJJx30sBN1/vnnFxw/8sgjatq0qc466yytXLmy4PohQ4ZoypQpkqS7775bLVq0UHJysu68805J0m+//aYBAwYoJSVFKSkp+uKLLyRJf//739WqVSu1atVKY8eOlSStXr1arVq1KvjeTz75pEaNGiXJjUQNHz5cHTp0UNOmTTVnzhzt27dPI0eO1KRJk9SmTRtNmjRJu3bt0tVXX60OHTqobdu2eu+99yRJe/bs0eWXX67mzZtrwIAB2rNnT5HPzYIFC9SlSxe1a9dOvXv31saNG5Wenq5mzZoV/NxXXHGFXnrpJUnS6NGj1b59eyUnJ+uBBx4o+DmaNWumQYMGqVWrVvrll1/0+OOPq3Xr1kpJSdHdd99dqudu8+bNuuiii9S+fXu1b9++xNLz9ddf64wzzlDbtm3VsWNHrVy5skzP0YQJE9S/f39169ZN3bt3lyRdcMEFmhiAf3Oy1oblR7t27SwAILosX7689Hd+801rq1SxVir8qFLFXV9BGRkZtmnTprZVq1Y2MzPzkNsbNGhgd+zYYa21Ni0tzbZq1cru2rXLpqen25NOOsmOHj3aWmvt4MGD7eTJk+3vv/9umzZtanNzc6211m7bts1aa+2ll15qx4wZY621Njs7227fvr3g+2VkZNidO3faFi1a2IULF9qff/7ZtmzZsiDD6NGj7QMPPGCttbZLly729ttvt9Za++GHH9ru3btba6199dVX7U033VTwNffcc4994403CjKcfPLJNiMjwz711FN26NCh1lprlyxZYuPj4+38+fMP+Jn37dtnzzjjDLtp0yZrrbVvvfVWwdd88skn9vTTT7f/93//Z3v37m2ttfbjjz+21157rc3NzbU5OTm2b9++dvbs2fbnn3+2xhg7b948a62106ZNs2eccYbdtWuXtdbaLVu2lOq5u+KKK+ycOXOstdauWbPGnnLKKYe8TjNnzrR9+/a11lqbnp5us7KyrLXWfvrpp/bCCy8s03P06quv2nr16hXks9badevW2VatWh3yuEX9O5aUZg/TadhnCwAQnkaMkHbvPvC63bvd9QMHVuhbV61aVZdddpmqVaumypUrH3L71q1bdeSRR0qS5syZowEDBqhKlSqSpP79+x9y/+rVqyspKUnXXHON+vXrV7COaMaMGXr99dclSfHx8apevbrmzp2rAQMGqGrVqpKkCy+8UHPmzCny++7vwgsvlCS1a9dOq1evLvI+n3zyiaZOnVqwpiwzM1Nr167V559/rltvvVWSlJycrOTk5EO+duXKlVq2bJl69uwpScrJyVHdunUlST179tTkyZN10003acmSJQWP9cknn6ht27aSpIyMDP3www9q0KCBGjZsqNNPP12SNH36dA0dOrTg+atVq1apnrvp06dr+fLlBffbsWOHMjIyVK1atSJ/9vT0dA0ePFg//PCDjDHKysoq03OU/3Pun+/YY4/Vhg0bivw+ZUHZAgCEp7z/AZb6+jKKi4tTXFzRq2kSEhKUm5t72NuLuv/XX3+tzz77TFOmTNGzzz6rGTNmlClP/mPmO3gvp/xSGB8fr+zs7CK/h7VW//73v9WsWbMyPXb+17Zs2VLz5s075Lbc3FytWLFCVapU0bZt21S/fn1Za3XPPffo+uuvP+C+q1evLiiSpXG45y43N1dffvmlkpKSSvV97r//fp199tl65513tHr1anXt2vWwP2dRz9FXX311SO7MzEwdccQRpf5ZDoc1WwCA8NSgQdmuD6BmzZrpp59+kiR17txZ7777rvbs2aOdO3fq/fffP+T+GRkZSk9PV58+fTRmzJiC0Z/u3btr/PjxktxIUXp6ujp16qR3331Xu3fv1q5du/TOO++oU6dOqlOnjjZt2qQtW7Zo7969+uCDD0rMeeSRR2rnzp0Fx71799a4cePkZrWkRYsWFfwM//rXvyRJy5Yt0zfffFPkz7x58+aCspWVlaVvv/1WkjRmzBg1b95c//rXvzR06FBlZWWpd+/eeuWVV5SRkSFJWr9+vTZt2nTI9+3Zs6deffVV7c4bpdy6dWupnrtevXpp3LhxBfdbvHhxsc9Fenq66tWrJ8mtvyrrc1SU77///oB1dOVF2QIAhKdHHpHypp4KVKnirg+yvn37atasWZKkU089VZdddplSUlJ07rnnqn379ofcf+fOnerXr5+Sk5N11lln6e9//7sk6emnn9bMmTPVunVrtWvXTsuXL9epp56qIUOGqEOHDjrttNP0xz/+UW3btlViYqJGjhypDh06qGfPnjrllFNKzHn22Wdr+fLlBYu/77//fmVlZSk5OVktW7bU/fffL0m64YYblJGRoebNm2vkyJFq167dId+rUqVKmjJlioYPH66UlBS1adNGX3zxhVauXKl//vOfeuqpp9SpUyd17txZDz/8sHr16qUrr7xSZ5xxhlq3bq2LL774gFKT75xzzlH//v2VmpqqNm3aFEzflfTcPfPMM0pLS1NycrJatGih559/vtjn4i9/+YvuuecetW3b9oCRv9I+R0WZOXOm+vbtW+zjlobJb3bhJjU11freAwQAEFgrVqxQ8+bNS/8FEye6NVpr17oRrUceqfB6rdLYuHGjBg0apE8//TToj4Xw1blzZ7333nuqWbPmAdcX9e/YGLPAWpta1PdhzRYAIHwNHBiScnWwunXr6tprr9WOHTtK3GsL0Wnz5s26/fbbDyla5UHZAgCgCJdeeqnvCPCodu3auuCCCwLyvVizBQAAEESULaAYXSd0VdcJXX3HAABEMMoWAABAEFG2gCLkj2jNXjNbs9fMZoQLQIHVq1cX7FklSWlpaQW7swfSu+++e8AO6iNHjtT06dMD/jgIPsoWAABlcHDZSk1N1TPPPBPwxzm4bD344IPq0aNHwB8HwUfZAoowa8gszRoyS10adlGXhl0KjgGE1sFbQQZia8jXX39dycnJSklJ0R/+8AetXr1a3bp1U3Jysrp3715wnrwhQ4bo1ltvVceOHXXiiSdqypQpkqS7775bc+bMUZs2bTRmzBjNmjWr4Hx+o0aN0tVXX62uXbvqxBNPLChhq1evPmAn8ieffFKjRo2SJL300ktq3769UlJSdNFFF2n37t364osvNHXqVN11111q06aNfvzxRw0ZMqQgw2effaa2bduqdevWuvrqq7V3715JUqNGjfTAAw/o1FNPVevWrfXdd99V/AlDhVG2AABhadQoadiwwoJlrTvO6yjl8u233+rhhx/WjBkztGTJEj399NO65ZZbNHjwYH3zzTcaOHDgAVOCGzdu1Ny5c/XBBx/o7rvvliQ99thj6tSpkxYvXqxhw4Yd8hjfffedPv74Y3399df661//etgTIue78MILNX/+fC1ZskTNmzfXyy+/rI4dO6p///4aPXq0Fi9erJNOOqng/pmZmRoyZIgmTZqkpUuXKjs7u+CUQJJ0zDHHaOHChbrhhhsO2a0dflC2gGIwogX4Ya20fbv09NOFhWvYMHe8fXv5R7hmzJihSy65RMccc4wkqVatWpo3b56uvPJKSdIf/vAHzZ07t+D+F1xwgeLi4tSiRQv99ttvpXqMvn37qnLlyjrmmGN07LHHlvh1y5YtU6dOndS6dWtNnDix4HyEh7Ny5Uo1btxYTZs2lSQNHjxYn3/+ecHtF154oSSpXbt2Wr16dakyI7jY1BQAEHaMkcaMcZefftp9SNKf/+yuNyY0OSpXrlxwubSnt9v/a+Lj45Wdna2EhATl5uYWXJ+ZmVlweciQIXr33XeVkpKiCRMmFJyTsaKZ8x8b/jGyBQAIS/sXrnwVLVrdunXT5MmTtWXLFknS1q1b1bFjR7311luSpIkTJ6pTp07Ffo8jjzyyyBMuF6dOnTratGmTtmzZor179+qDDz4ouG3nzp2qW7eusrKyNHHixBIfp1mzZlq9erVWrVolSXrjjTfUpUuXMuVBaFG2AABhKX/qcH/7r+Eqj5YtW2rEiBHq0qWLUlJSdPvtt2vcuHF69dVXlZycrDfeeENP5w+jHUZycrLi4+OVkpKiMQe3wcNITEzUyJEj1aFDB/Xs2VOnnHJKwW0PPfSQTjvtNJ155pkHXH/55Zdr9OjRatu2rX788ceC65OSkvTqq6/qkksuUevWrRUXF6c//elPZXwmEEqmtMOioZaammrT0tJ8xwAABNCKFSvUvHnzEu+3/xqt/KnDg49DNZUIHKyof8fGmAXW2tSi7s+aLQBA2DFGqlHjwGKVP4hUowZFC5GFsgUACEujRrkRrvxilV+4KFqINKzZAgCErYOLFUULkYiyBQAAEESULQAAgCCibAEAAAQRZQsAgCDp06ePtm/f7jsGPOPdiAAAVED+6XiKMm3atBCnQThiZAsAEFPefPNNdejQQW3atNH111+vnJwczZ8/X8nJycrMzNSuXbvUsmVLLVu2TLt27dLVV1+tDh06qG3btnrvvfckSRMmTFD//v3VrVs3de/eXRkZGRo6dKhat26t5ORk/fvf/5YkNWrUSL///rt27dqlvn37KiUlRa1atdKkSZMkSQsWLFCXLl3Url079e7dWxs3bjwk75AhQ3TrrbeqY8eOOvHEEzVlyhRJ7lyNd911l1q1aqXWrVsXfM+bbrpJU6dOlSQNGDBAV199tSTplVde0YgRI4L75KJIjGwBAMJa1wldJUmzhsyq8PdasWKFJk2apP/9739KTEzUjTfeqIkTJ2rQoEHq37+/7rvvPu3Zs0dXXXWVWrVqpXvvvVfdunXTK6+8ou3bt6tDhw7q0aOHJGnhwoX65ptvVKtWLQ0fPlzVq1fX0qVLJUnbtm074HE/+ugjHX/88frwww8lSenp6crKytItt9yi9957T7Vr19akSZM0YsQIvfLKK4fk3rhxo+bOnavvvvtO/fv318UXX6z//Oc/Wrx4sZYsWaLff/9d7du3V+fOndWpUyfNmTNH/fv31/r16wsK3Jw5c3T55ZdX+DlE2VG2AAAx47PPPtOCBQvUvn17SdKePXt07LHHSpJGjhyp9u3bKykpSc8884wk6ZNPPtHUqVP15JNPSpIyMzO1du1aSVLPnj1Vq1YtSdL06dMLTmYtSTVr1jzgcVu3bq077rhDw4cPV79+/dSpUyctW7ZMy5YtU8+ePSVJOTk5qlu3bpG5L7jgAsXFxalFixb67bffJElz587VFVdcofj4eNWpU0ddunTR/Pnz1alTJ40dO1bLly9XixYttG3bNm3cuFHz5s0r+LkQWpQtAEBYyh/Rmr1m9gHHFRnhstZq8ODBevTRRw+5bcuWLcrIyFBWVpYyMzNVtWpVWWv173//W82aNTvgvl999ZWqVq1a6sdt2rSpFi5cqGnTpum+++5T9+7dNWDAALVs2VLz5s0r8esrV658wM9QnHr16mn79u366KOP1LlzZ23dulVvv/22qlWrpiOPPLLUmRE4rNkCAMSM7t27a8qUKdq0aZMkaevWrVqzZo0k6frrr9dDDz2kgQMHavjw4ZKk3r17a9y4cQUFZ9GiRUV+3549e+q5554rOD54GnHDhg2qUqWKrrrqKt11111auHChmjVrps2bNxeUraysLH377bel/lk6deqkSZMmKScnR5s3b9bnn3+uDh06SJJOP/10jR07tmBa8cknn1SnTp1K/b0RWIxsAQDCUv4IViDXbLVo0UIPP/ywevXqpdzcXCUmJuq5557T7NmzlZiYqCuvvFI5OTnq2LGjZsyYofvvv1+33XabkpOTlZubq8aNG+uDDz445Pved999uummm9SqVSvFx8frgQce0IUXXlhw+9KlS3XXXXcpLi5OiYmJGj9+vCpVqqQpU6bo1ltvVXp6urKzs3XbbbepZcuWpfpZBgwYoHnz5iklJUXGGD3xxBM67rjjJLki9sknn6hJkyZq2LChtm7dStnyyJQ0HOlLamqqTUtL8x0DABBAK1asUPPmzcv0NYEsW0AgFPXv2BizwFqbWtT9GdkCAIQ1ShYiHWu2AAAAgoiyBQAAEESULQBASIXrWmGgNMrz75eyBQAImaSkJG3ZsoXChYhkrdWWLVuUlJRUpq9jgTwAIGTq16+vdevWafPmzb6jAOWSlJSk+vXrl+lrKFsAgJBJTExU48aNfccAQoppRAAAgCCibAEAAAQRZQuIEF0ndC3YSRsAEDkoWwAAAEHEAnkgzOWPZs1eM/uAY05hAgCRgZEtAACAIGJkCwhz+SNYjGgBQGRiZAsAACCIGNkCIgQjWgAQmRjZAgAACCLKFgAAQBBRtgAAAIKIsgUAABBElC0AAIAgomwBAAAEEWULAAAgiChbAAAAQUTZAgAACCLKFgAAQBBRtgAAAIKIsgUAABBElC0AAIAgomwBAAAEEWULAAAgiChbAAAAQUTZAgAACCLKFgAAQBBRtgAAAIKIsgUAABBElC0AAIAgomwBAAAEUUDKljHmHGPMSmPMKmPM3UXcfrsxZrkx5htjzGfGmIaBeFwAAIBwV+GyZYyJl/ScpHMltZB0hTGmxUF3WyQp1VqbLGmKpCcq+rgAAACRIBAjWx0krbLW/mSt3SfpLUnn738Ha+1Ma+3uvMMvJdUPwOMCAACEvUCUrXqSftnveF3edYdzjaT/BuBxAQAAwl5CKB/MGHOVpFRJXQ5z+3WSrpOkBg0ahDAZAABAcARiZGu9pBP2O66fd90BjDE9JI2Q1N9au7eob2StfdFam2qtTa1du3YAogEAAPgViLI1X9LJxpjGxphKki6XNHX/Oxhj2kp6Qa5obQrAYwIAAESECpcta222pJslfSxphaS3rbXfGmMeNMb0z7vbaEnVJE02xiw2xkw9zLcDwsfEiVKjRlJcnPs8caLvRACACBSQNVvW2mmSph103cj9LvcIxOMAITNxonTdddLuvDfRrlnjjiVp4EB/uQAAEYcd5IGijBhRWLTy7d7trgcAoAwoW0BR1q4t2/UAABwGZQsoyuG2HmFLEgBAGVG2gKI88ohUpcqB11Wp4q4HAKAMKFtAUQYOlF58UWrYUDLGfX7xRRbHAwDKLKQ7yAMRZeBAyhUAoMIY2QIAAAgiyhYAAEAQUbYARA529QcQgVizBSAysKs/gAjFyBaAyMCu/gAiFGULQGRgV38AEYqyBSAysKs/gAhF2QIigLXFH8cEdvUHEKEoW0CYGzVKGjassGBZ645HjfKZygN29QcQoShbQBizVtq+XXr66cLCNWyYO96+PQZHuAYOlFavlnJz3WeKFoAIQNkCDqOiU3cZGe5Dkvbtk95+W1qxwh3v2uXeRPfFF+44PV268Ubp88/d8fbt0r33SgsXSmPGSDfc4ApWXJz7fMst0n33lftHAwCEEGULKEJJU3fWSqtWua2eJCknR7rySum119xxZqZ05JGuGOXfftll0tSp7njfPumJJ6QFC9zx3r3S5MlusEZy5Wv0aGnpUjdjdvPNB+b74x+l2rXd10jS999L559f+P127nRfu3dvAJ8UAEC5sKkpcJD9p+4kN7KUP3XXs2dhATvtNOmSS6Tnn5fi4135Ou00d1tSkvTUU9KZZxYeL10q1avnjmvWlLKyCh/z2GOlzZsLjxs2LLzdWmn8+AMzjhvncqWmuuOMDOmnn9zsmiTNnSv16eM+n3mmtGSJ9Oab0u23S3XrBuRpAgCUkrFhuugjNTXVpqWl+Y6BGJJfsmrWdJdbtZKWLy+8/bjjpB49pDfecMcffiidcIKUnBzcTPlF789/PrD45R8bc+jX/fqrNHu2dM45UvXqbvP1q6+Wfv5ZOv546fXXpbFjpY8/diNkO3a4N/Yl8OcXAJSLMWaBtTa1qNuYRkTMys11o0H5Bg2SOnZ0l42RLr30wPv/8kth0ZKkvn2DW7Tyc9SocWCxGjPGHdeoUXTRklwxvOwyV7Qkt458167CUa3q1aX69aWjj3bHf/ubVKdO4Wjahg3Snj3B/MkAIHbwdyxiyoYNrnAYI919t/TMM259VOXKbkqwY8fCacJt2w782jvvPPxIUjCNGuUy5T9ufuEqa479R63OP9995OvdW6pVS0pMdMc33yx9+620cqU7/vVXN9UZx59nAFBmlC1Etexs9zkhwY1KDRok/fCD1KSJW9DeqlVhuerf330ubupO8lO4Dn68QD/+2We7j3w33uimVPN16+aeq7ffdsd797qCCgAoGWULUWvpUqlrV7c+qW9fqUsXt2j9qKPc7W3auI+DHW7qTip+6i6a9OhReNlaNwpYu7Y73rPHLfT/61/dFhQAgOKxQB5RIzNTGjLEjdBcf70bfbnhBjdKk1rkksXi7T91V9RxrNq6VXr0UTcNedZZbt3boEFuSvbUU32nAwA/ilsgz8gWItrMmW490RVXuO0Vtm51C8ElN831yivl/97BnrqLVLVquT3A8m3a5Na91arljhculJYtc28wSErykxEAwgkjW4g4mzcXTmlddJH7H/t331GGwsWdd0ovvOBK2BFHSFu2uCLG6wMgmrH1A6LGP/7htizYtMkdP/us27CT/5GHj9GjpcWLXdGSXCHu29drJADwirKFsLZli3TPPW70SnLvirv//sJtDOrWZaoq3BgjnXSSu2yt20z1qqvccW6ue+PBokXl//5dJ3RV1wldK5wTAEKFNVsIO/vv5C650auGDd3WA6ecwgmYI4kxbvF8vlWr3PkjTz9datvWvYkhJ8ftXg8A0YqyhbDTo4f7n+/777sdztetK9wJHZGtaVNp/XqpUiV3/MYb0vDh7gTajRoV/7X5o1mz18w+4HjWkFlByQoAgcI0IrxLT5defbVwc9ErrjjwVDkUrehStWrhTvUpKW67joYN3fFHH7nRLwCIJrwbEd49/7zbD2vhQje1hNiUmys1biy1bClNm3b4+zGiBSAc8W5EhJVdu6S77nLThJI0eLCUlkbRinVxcdJXX0ljx7rjrVulyy8vPD+jVDj6ebhjAAhHrNlCyOTkSPHx7t2DH34oVasmnXee2yKgXTvf6RAOjjvOfUjSN99In34qjRjhjkeOlHbskGaOmSVjCs9hWaOGO1k3AIQryhZCYvx4t9Hl/Pluvc7ChWzZgOJ17ereHHHEEa5Yvf++27/LWjf6tf/JwjmVEoBwRtlC0Oze7UayKleWGjRwa3F27nS7iVO0UBr5G6Ma4944YYw7B+Mzz7jr9z9ZOACEK9ZsISh+/VVq0sQtfpfcDuITJxaePw8oq7/8xW0Rsb9OnShaAMIfZQsBtX69+3zccdKVV0rt2/vNg+iRv0Zrfx995K7//Xdp3z4/uQCgJJQtBMxf/yq1bu3eRSZJTz4pdezoNxOiQ37Ryl+jlX/an3/+011/zTXu31puru+kAHAo1myhQrZtc+uyjjrKnXC4alX3AQSSMe5dh/uv0Rozxt1Wo4Z0/fXSb7+57SMk6ZdfpBNO8JUWAA7EpqYot5073elXLrmkcMEyEEwHv+uwqHchfvaZdM45bnuRXr1Cmw9A7GJTUwTUpk3u85FHSvfcIw0d6jcPYsfBxaqoxfFt2kh33il17uyOf/+dzU8B+EXZQpm89ZY7j92KFe741lvZ+R3h5eijpUcfdduLZGe7E5vzBwEAn1izhVLZu9ftl9W9u/SnP0l16vhOBJTMGLeeq149d2ytK2D5J8IGgFBgZAsluuYat/jdWql2bbcwmf2yEAni491Jzvv3d8cvvyydeqpbTA8AocLIFkrUtq1712FurvufFxCp6tWTUlKkY4/1nQRALGFkC4f4/XepTx9p+nR3fPPN0v33U7QQ+c49V3rzTTe9mJ4unXWW9MUXvlMBiHaULRyialU3zbJxo+8kQPBs2CBt2cL6LQDBR9mCJHcuw7vukrKy3Ml/58+X/vAH36mA4GneXFq2rPCUUuPHM8oFIDgoW5Dk/ifz7LOFJ/qN418GYkD+1HhmpvT3v0svvOA3D4DoxP9SY1hGhjRvnrt84YXSjz9Kp5/uNxMOY+JEqVEj14IbNXLHCJikJPeHRv6ZENavl777zm8mANGDshXDrr9e6tfPnXZHko4/3m8eHMbEidJ110lr1rj9N9ascccUroA66iipenV3+bbb3A70u3d7jQQgSnBuxBiTk+PWZSUlST/9JK1bV3haE4SpRo1cwTpYw4bS6tWhThMTNmyQlixx716UXOmqUsVvJgDhjXMjQpLbObtXL3eKHUk68USKVkRYu7Zs16PCjj++sGh9+KF08smFp6gCgLKibMWQhASpa1fpzDN9J0GZNGhQtusRUPXquf24Gjf2nQRApKJsRbnsbGnkSOmbb9zx/fdLgwf7zYQyeuSRQ+ewqlRx1yPo2rSRJk1yU+9797r1XJs2+U4FIJJQtqLctm3SSy9J777rOwnKbeBA6cUX3RotY9znF1901yOk0tLcU//1176TAIgkLJCPUt9+K7Vo4f7fvGkT54IDAuXXX6XjjnOXV6yQTjnF/XcGILaxQD7GfPmlO9nuhAnumKIFBE5+0Vq71u0+P2qU1zgAIgBlKwp16CA99ph08cW+kwDR64QTpMcfd/vVlQsb1QIxg7IVJVaulPr0cSfWjYuT7rxTOvJI36mA6GWMdNNNbpsIa6Ubb5QmTy7lF7NRLRBTKFtRIj3dnVT35599JwFiz65d0qJFZTjFz4gRh25Pv3u3ux5A1EnwHQDlZ617d1T79m7qcNUqqVIl36mA2FOtmjR7ttvLTpJ++EGqU8edAqhIbFQLxBRGtiLYc8+5E0cvWuSOKVqAP5UquSn8rCw3pV/smkk2qgViCiNbEWzIECkx0W26CCA8JCZKL7zgRrsO65FH3Bqt/acS2agWiFqMbEWYBQukq65yfz1Xq+beCcUeP0B46dbNTe1L0tix0vPPH3QHNqoFYgplK8IsXy7NnStt2OA7CYCSWCvNnClNn+4uH2DgQGn1aik3132maAFRix3kI4C10i+/FC7n2L370FPlAQhP2dnuIynJnT6rUiWpalXfqQAEGjvIR7jHHnM7wq9Z444pWkDkSEhwRctat2i+d283mHXw37lh+ncvgABggXwEuOIK95fxCSf4TgKgvIyRhg1ze3I9+KC0fbs0Zoy73lp3W40anP4HiEaMbIWpX391v4itdWfyuP9+97ZyAJGrXz/p0ktd0Xr6aTfSlV+0nn7aXc8IFxB9WLMVph56yJ13belSqXFj32kABFJurlSvnvujKt+f/1w40gUg8hS3ZouyFWZyc90IVm6u2xG+aVPfiQAEw9at0tFHFx7n5lK0gEjGAvkI8f777tQ7W7e6wkXRAqKTtW7d1v7atZNycg69H4DIR9kKI9Wru41KD/6FCyB67L9G689/dm9+SUlxp91KTXUjXPvfjwXzQOSjbIWBFSvc586dpVmzpNq1vcYBEETGuHcd5q/Rio93Z4ZISZEWL5auvda9Y5FF80D0YM2WZy+/LP3pT9K8ee6vWgCxwdoD12jl5kq33upOMJ+PRfNA5GCBfBjbsUMaP1668073Fy6A2NF1QldJ0qwhsyS5Arb/Fi8smgciBwvkw8y+fdJTT7m1GkcdJQ0fTtECYl3+Gq399ekj/fyznzwAAocd5D2YNs2NZLVsKZ1zju80AEItf0Rr9prZBcerVknrn55VMHV4443S889LfftK337LCBcQyShbHlxwgVsIm5LiOwmAcJGQcOAarX/8wy2UP+EEihYQ6VizFSLp6dLQodLf/iadcorvNADCQVFrtvYvVvnHe/e63x933imdemrocwIoGWu2wsDmzdL8+dLKlb6TAAhXxrgCll/C8ovXpk3uHctLl/rLBqD8mEYMssxMKSlJatJE+v576YgjfCcCEC7yR7RKcsIJbt1WlSruOP/3CoDIwMhWEG3Z4vbOeuYZd0zRAnA4+SNas9fM1uw1sw8Y4ZIKi9Y337g/3j7/3E9OAGVH2Qqio45y5ztLTvadBEC0OO449+aaBg18JwFQWiyQD4JNm9wo1pFH+k4CINIcvGi+JD//LDVuHLw8AEqHBfIhlJ0t9ewpXXIJ5zMDEFz/+Ifbr+/bb30nAVAcFsgHWEKCNHKkdMwx7I0DoOxKO6IlSRdf7N7p3Lx58PIAqDimEQNk82bpxx+l00/3nQRALNq6VfrhB+m003wnAWIT04ghcPPN0nnnSRkZvpMAiEU33CD168fvICAcBaRsGWPOMcasNMasMsbcXcTtlY0xk/Ju/8oY0ygQjxtOnnlGmjJFqlbNdxIAsWjMGOntt/kdBISjCpctY0y8pOcknSuphaQrjDEtDrrbNZK2WWubSBoj6fGKPm44yMiQxo1zC+Hr1JG6dPGdCEBEmzhRatRIiotznydOLPWXHn+8dPbZ7vLHH0tLlgQlIYByCMTIVgdJq6y1P1lr90l6S9L5B93nfEmv5V2eIqm7MZG/fPz116Vhw6RFi3wnARDxJk6UrrtOWrPG/QW3Zo07LkPhkqR9+6SbbpKGDw9STgBlFoiyVU/SL/sdr8u7rsj7WGuzJaVLOjoAj+3VDTdIaWmcGBZAAIwYIe3efeB1u3e768ugUiU3svXWWwHMBqBCwmqBvDHmOmNMmjEmbfPmzb7jFCk7W7rrLmn9ere1Q5s2vhMBiApr15bt+mKcdJJUo4aUkyPde6/7fQXAn0CUrfWSTtjvuH7edUXexxiTIKm6pC0HfyNr7YvW2lRrbWrt2rUDEC3wvv9eeuEFado030kARJXDnX+nAuflWbVKevZZ6b33yv0tAARAIDY1nS/pZGNMY7lSdbmkKw+6z1RJgyXNk3SxpBk2XDf4KkGLFtJ337nFqAAQMI884tZo7T+VWKWKu76cmjWTVqyQ6h28sANASFV4ZCtvDdbNkj6WtELS29bab40xDxpj+ufd7WVJRxtjVkm6XdIh20OEu7/9Tfq//3OXKVoAAm7gQOnFF6WGDd0ahYYN3fHAgRX6tvlFa9UqadAgac+eAGQFUCbsIF8KWVlSr17uZK+vvOI7DQCU3eTJbvPlWbM4vQ8QDMXtIM+5EUshMdGt0UpM9J0EAMrnkkuk3r2lo47ynQSIPWH1bsRw8+mn0mWXuWH3I45wJ5kGgEiVX7TGj5fuucdvFiCWUB+K8eOP7t2H+/a5sgUA0WDZMrejRHY2f0QCocCarSJY69anSq5oVarkJQYABEVOjvs9l5Bw4O87AOVX3JotphEPsnGjdMYZ0oIF7piiBSDaxMe7orVjh3TuuW7HeQDBQ9k6SEaGtGuX7xQAEHzWSlu3Sr//7jsJEN2Yrc+TP5R+8snSkiVSHDUUQJSrXl2aN8+NdElMKQLBQqWQ+wUzdKg0apQ7pmgBiBX5RWv2bLeEYsshJ1IDUFHUCrnFonFx/EUHIHbFx7sNnPc/WxCAwIj5acTcXLdQ9OWXfScBAH/OOkuaP9/94Zn/JnX+AAUCI6ZHtv77Xzds/ttv7pcKv1gAxLK4OPcH6J//LN1xh+80QPSI6bJlrVS5slSliu8kABAe9v/DM0y3YQQiTsxvasq7bwDgQPv/XuR3JFA6bGpaDH6JAMCB8n8v/vyzdNppbjscAOUX82ULAFC0ypWlvXvdTvMAyi/m340IACja8cdLixYV7j3IlCJQPoxsAQAOK79ovfqqdPnl7t2KAMqGsgUAKNGOHdL27dKePb6TAJGHsgUAKNGtt0rTpklVq/pOAkQeyhYAoETGuFP67NghXXSRO4E1gNKhbAEASi0rS1q+XFq50ncSIHLwbkQAQKkdfbTbd6tSJd9JgMjByBYAoEzyi9bcudL11/MORaAklC0AQLl8/bU0c6a0ebPvJEB4o2wBAMpl2DBp4UKpTh3fSYDwRtkCAJSLMVK1alJOjvTAA9K33/pOBIQnFsgDACrk99+lF15wu823bOk7DRB+KFsAgAqpU0davFg67jjfSYDwxDQiAKDC8ovW2rXS+PF+swDhhrIFAAiYceOke+6RfvvNdxIgfFC2AAAB88gj0oIFvEMR2B9lCwAQMJUqSSed5C6/+66Unu41DhAWKFsAgIBbvVq65BJp9GjfSQD/eDciACDgGjWSPv1U6tjRdxLAP0a2AABB0bWrm1bcvdttDQHEKsoWACCorrlG6tVLysjwnQTwg7IFAAiqBx+UXn/dndoHiEWs2QIABNXJJ7sPSVq/XqpXz28eINQY2QIAhMRXX0lNmkhvv+07CRBalC0AQEiceqp0yy1u4TwQS5hGBACERGKi9MQT7rK10r59UuXKfjMBocDIFgAgpKyVLrvMvUvRWt9pgOBjZAsAEFLGSG3auD24gFhA2QIAhNy99/pOAIQO04gAAG++/FLq0UPaudN3EiB4KFsAAG/27ZPWrpU2bPCdBAgephEBAN507iwtXy4l8H8jRDFGtgAAXiUkSDk50t/+Ji1b5jsNEHiULQCAd1u3Ss88I731lu8kQOAxcAsA8K52bWnhQqluXd9JgMBjZAsAEBaOP97twbV+vTR9uu80QOAwsgUACCs33ijNny/99JOUlOQ7DVBxlC0AQFgZN07as4eihehB2QIAhJUGDQovf/+91LSpvyxAILBmCwAQliZNkpo3l+bO9Z0EqBjKFgAgLPXrJ/31r1Jqqu8kQMUwjQgACEtVq0r33ecu79snJSa6dysCkYaRLQBAWNu6VTr9dOm553wnAcqHsgUACGs1a0qtWkkNG/pOApQP04gAgLBmjPT6675TAOXHyBYAIGK89ppbNA9EEka2AAAR48svpZUrpexsKYH/gyFC8E8VABAxxoxx70qMj/edBCg9phEBABEjKckVrR07pH/8Q7LWdyKgZJQtAEDEef116ZZbpCVLfCcBSkbZAgBEnBtukBYskNq08Z0EKBllCwAQceLjC4vWihVSbq7XOECxKFsAgIj1zTdSSoo0frzvJMDhUbYAABGrdWvp0Uelyy/3nQQ4PLZ+AABELGOkO+5wl62VcnLYfwvhh5EtAEDE27tX6tNHGjXKdxLgUJQtAEDEq1xZOvFEqV4930mAQzHYCgCICs895zsBUDRGtgAAUeXjj6W//c13CqAQZQsAEFXef1/617+kPXt8JwEcphEBAFHliSekuDh3HkUgHDCyBQCIKlWquKK1b5/0wQe+0wCULQBAlHrmGem886SlS30nQaxjGhEAEJVuuklq2dLtMg/4xMgWACAqHXGEdO657vL27V6jIMZRtgAAUe3rr6VGjdyWEIAPlC0AQFRLTpYuukhq0sR3EsQq1mwBAKJaUpL08su+UyCWMbIFAIgJO3dK114rzZjhOwliDWULABATEhKk//1PWrzYdxLEGqYRAQAx4YgjpIUL2VkeocfIFgAgZuQXrcWL3bsUgVBgZAsAEFNycqTLL5dq15bmzPGdBrGAsgUAiCnx8dLkyVK9er6TIFYwjQgAiDmtW0u1aknWSlu2+E6DaEfZAgDErKFDpe7dpaws30kQzZhGBADErIsuktq3d1OLQLBUqGwZY2pJmiSpkaTVki611m476D5tJI2XdJSkHEmPWGsnVeRxAQAIhPPO850AsaCi04h3S/rMWnuypM/yjg+2W9Iga21LSedIGmuMqVHBxwUAIGA++US68kopN9d3EkSjipat8yW9lnf5NUkXHHwHa+331tof8i5vkLRJUu0KPi4AAAGzYYPbe+vXX30nQTSqaNmqY63dmHf5V0l1iruzMaaDpEqSfqzg4wIAEDCDB7uydfzxvpMgGpVYtowx040xy4r4OH//+1lrrSRbzPepK+kNSUOttUUO1BpjrjPGpBlj0jZv3lzGHwUAgPIxRqpUSdq3T3rrLbclBBAoJS6Qt9b2ONxtxpjfjDF1rbUb88rUpsPc7yhJH0oaYa39spjHelHSi5KUmprKP3UAQEhNmCBdf73UsKF0xhm+0yBaVHQacaqkwXmXB0t67+A7GGMqSXpH0uvW2ikVfDwAAILm6qul6dMpWgisipatxyT1NMb8IKlH3rGMManGmH/m3edSSZ0lDTHGLM77aFPBxwUAIOASEtwmp5K0c6ffLIgexobpxHRqaqpNS0vzHQMAEIPmznV7cH34odSxo+80iATGmAXW2tSibuN0PQAAHKRNG6lfP6lOse+xB0qH0/UAAHCQatWkN97wnQLRgpEtAAAOY9s26U9/kn74wXcSRDLKFgAAh5GZKU2ZIs2Z4zsJIhnTiAAAHEbdutJPP0lHHeU7CSIZI1sAABQjv2gtWeKmFYGyomwBAFCCX3+VOnSQHn7YdxJEIqYRAQAowXHHSW++WbjhKVAWlC0AAErhkkvcZ2ul3FwpPt5vHkQOphEBACilXbukXr2ksWN9J0EkoWwBAFBKVaq4XeWrV/edBJGEaUQAAErJGLd2CygLRrYAACiHDz+UZs70nQKRgJEtAADKKDtbuuMOqUkT6eyzfadBuKNsAQBQRgkJ0rRpUv36vpMgEjCNCABAOZx4olSpkrRvn7R5s+80CGeMbAEAUE7WSp07SzVrupEuY3wnQjiibAEAUE7GSDfeKNWqRdHC4VG2AACogEGDfCdAuGPNFgAAAfDcc9Kjj/pOgXDEyBYAAAEwf75bKJ+bK8UxlIH9ULYAAAiA55+XKldm7RYORfcGACAAkpJc0dq6VVq40HcahBNGtgAACKCLL5bWrpVWrpTi432nQTigbAEAEEBPPOE2O6VoIR9lCwCAAEpNLbzMYnlIrNkCACAoHnhAuuQS3ykQDihbAAAEwVFHSUcfLWVl+U4C35hGBAAgCO64w3cChAtGtgAACKIffpDeecd3CvhE2QIAIIjuvlu6+WZp3z7fSeAL04gAAATR2LFSYqLbDgKxibIFAEAQnXBC4eWMDKlaNX9Z4AfTiAAAhMANN0jdurm9txBbGNkCACAEzj5batBAyslho9NYQ9kCACAELr3UdwL4QrcGACCEZs2Sxo3znQKhRNkCACCE3nzTla29e30nQagwjQgAQAg9+aRUubL7QGxgZAsAgBCqUUM64gi3UH79et9pEAqULQAAPLjsMql3byk723cSBBvTiAAAePDHP0rbt0vx8b6TINgoWwAAeHDOOb4TIFSYRgQAwKPXXpMefth3CgQTZQsAAI+++EL65BO3YB7RiWlEAAA8GjtWSkqSjPGdBMHCyBYAAB4dcYQrWhkZ0vff+06DYGBkCwCAMHDOOdLOndKiRZyoOtpQtgAACAMPPuhGuSha0YeyBQBAGOjWzXcCBAv9GQCAMGGtNHIkW0FEG0a2AAAIE8ZIP/8sJSa64sU7FKMDZQsAgDDy6qtSAv93jipMIwIAEEbyi9Yvv0jffec3CwKD7gwAQJjJyZG6dJEaN5Y++8x3GlQUZQsAgDATHy+9/LJ04om+kyAQKFsAAIShs8/2nQCBwpotAADCVGamdNVV0rPP+k6CiqBsAQAQpipXltLTpd27fSdBRTCNCABAmDJGmjqV/bYiHSNbAACEsfyiNX++2w4CkYeyBQBAmNu2zW0F8cgjvpOgPJhGBAAgzNWsKb33nnT66b6ToDwoWwAARICePd1na91n1nFFDqYRAQCIEBs2SJ06SdOm+U6CsqBsAQAQIWrXluLi3P5biBxMIwIAECESE6XPP/edAmXFyBYAABHGWumdd6Q9e3wnQWlQtgAAiDBpadKFF0qvveY7CUqDaUQAACJM+/bSRx9JPXr4ToLSoGwBABCBevd2n61lG4hwxzQiAAAR6ssvpRYtpNWrfSdBcShbAABEqPr13e7y27b5ToLiMI0IAECEql9f+uIL3ylQEka2AACIcJmZ0uTJvlPgcChbAABEuH/+U7r0UmnxYt9JUBSmEQEAiHB//KPUurXUpo3vJCgKI1sAAES4pCSpSxd32Vq/WXAoyhYAAFHilVekjh2l7GzfSbA/yhYAAFGiZk3p2GOl9HTfSbA/1mwBABAlBgxwHwgvjGwBABBlfvtN+uwz3ymQj5EtAACizA03SP/7n7R2rVS5su80oGwBABBlHn/cvSuRohUeKFsAAESZk08uvGytZIy/LGDNFgAAUcla6cYbpTvv9J0ElC0AAKKQMVJcnPuAX0wjAgAQpcaNYwoxHNB3AQCIUvlFa9ky6aef/GaJZZQtAACi2K5d0llnSaNG+U4Su5hGBAAgilWtKk2eLLVr5ztJ7KJsAQAQ5Xr29J0gtlVoGtEYU8sY86kx5oe8zzWLue9Rxph1xphnK/KYAACg7Navl3r1kmbP9p0k9lR0zdbdkj6z1p4s6bO848N5SNLnFXw8AABQDrVqSZs2uQ+EVkWnEc+X1DXv8muSZkkafvCdjDHtJNWR9JGk1Ao+JgAAKKMjjpAWLWIrCB8qOrJVx1q7Me/yr3KF6gDGmDhJT0liD1sAADwyxu0s//nn7jNCo8SyZYyZboxZVsTH+fvfz1prJRX10t0oaZq1dl0pHus6Y0yaMSZt8+bNpf4hAABA6bz/vtSli/Thh76TxI4SpxGttT0Od5sx5jdjTF1r7UZjTF1JRc0EnyGpkzHmRknVJFUyxmRYaw9Z32WtfVHSi5KUmppK5wYAIMD69JEmTHCL5REaFV2zNVXSYEmP5X1+7+A7WGsH5l82xgyRlFpU0QIAAMGXkCANHuw7RWyp6JqtxyT1NMb8IKlH3rGMManGmH9WNBwAAAiO6dOlCy6QsrN9J4l+FRrZstZukdS9iOvTJP2xiOsnSJpQkccEAAAVl5EhrVwprVsnNWrkO01049yIAADEoPPPdyeopmgFH2ULAIAYZIwUHy9lZUkrVvhOE90oWwAAxLBBg6QePaTMTN9JohcnogYAIIbddpt01VVS5cq+k0QvyhYAADHstNN8J4h+TCMCABDjcnOlMWOkf7JpU1BQtgAAiHFxcdJ//yvNmOE7SXRiGhEAAOidd6SqVX2niE6MbAEAgIKitW2b2/AUgUPZAgAAkqRNm6QTT5TGjvWdJLowjQgAACRJxx4r/eUvUt++vpNEF8oWAAAocM89vhNEH6YRAQDAATZtciNc27f7ThIdKFsAAOAAGze6fbdmzvSdJDowjQgAAA6QkiL98ot03HG+k0QHRrYAAMAh8osW20BUHGULAAAU6YUXpEaN3N5bKD/KFgAAKNIZZ0gDB7pzJ6L8WLMFAACKlJwsPf207xSRj5EtAABQrKVLpbff9p0iclG2AABAsf76V+mOO6SsLN9JIhPTiAAAoFhjxrgTVScm+k4SmShbAACgWCecUHg5N1eKY16sTHi6AABAifbtk3r3lh56yHeSyEPZAgAAJapUSWrYUKpd23eSyMM0IgAAKJUXX/SdIDIxsgUAAErNWunjj6WdO30niRyULQAAUGrLlknnnCO9/LLvJJGDaUQAAFBqrVtLH3wg9ezpO0nkoGwBAIAy6dvXd4LIwjQiAAAos1mzpG7dpF27fCcJf5QtAABQZomJ0q+/SmvX+k4S/phGBAAAZXbmmW6xPLvJl4ynCAAAlEtcnDs59YoVvpOEN8oWAAAotz/8QerVy53OB0VjGhEAAJTbbbdJgwa5NVwoGmULAACU2+mn+04Q/phGBAAAFZKVJT31lPSf//hOEp4oWwAAoEISEqQ33pA++sh3kvDENCIAAKgQY6TZs6Xq1X0nCU+MbAEAgArLL1rp6VJurt8s4YayBQAAAmLJEqlBA+n9930nCS+ULQAAEBAtW0pXXSU1aeI7SXhhzRYAAAiIhATpued8pwg/jGwBAICA2rBBGjfOd4rwQdkCAAABNXmyNGyY9P33vpOEB8oWAAAIqGuvdUWraVPfScIDZQsAAARUlSrSiSe6y9b6zRIOKFsAACAohg2TrrzSdwr/eDciAAAIimOOcZ+tdbvMxyrKFgAACIoRI3wnCA9MIwIAgKBaulRat853Cn8oWwAAIGi2bZPat5cef9x3En+YRgQAAEFTs6bbd+vMM30n8YeyBQAAguq883wn8ItpRAAAEHRLl0oDBrhpxVhD2QIAAEGXmyvNmyctX+47SegxjQgAAIIuJUX65RcpMdF3ktBjZAsAAIREYqLb4HTDBt9JQouyBQAAQuaGG6QzzpCysnwnCR2mEQEAQMhccYXUrp3vFKFF2QIAACHTpYv7iCVMIwIAgJDKyZHeekuaMcN3ktCgbAEAgJCy1p2k+qWXfCcJDaYRAQBASCUkuFGtE07wnSQ0GNkCAAAh17ChFBfnphSjHWULAAB4sWCB1KSJtGiR7yTBRdkCAABeNGkiNWsmZWf7ThJcrNkCAABeVK8uffSR7xTBx8gWAADwKiMjureBoGwBAACv7r1X6tNH2rrVd5LgoGwBAACvbr9dmjVLqlXLd5LgYM0WAADwqlEj9xGtGNkCAADe5eS4XeWff953ksBjZAsAAHgXHy999ZWUnu47SeBRtgAAQFj473+lxETfKQKPaUQAABAW8ovWb79Jubl+swQSZQsAAISNr76SGjSQPvzQd5LAoWwBAICwceqp0rBhUsuWvpMEDmu2AABA2EhMlB57zHeKwGJkCwAAhJ2VK6Xx432nCAzKFgAACDtvvCHdeae0ZYvvJBVH2QIAAGHnjjukn3+Wjj7ad5KKY80WAAAIOzVrFl62VjLGX5aKYmQLAACEpexsacAA6YEHfCepGMoWAAAISwkJ0jHHSDVq+E5SMUwjAgCAsPXSS74TVBwjWwAAIKxZK33xhZtWjESULQAAENZmzpTOPFP6z398JykfyhYAAAhrXbtKr74q9evnO0n5sGYLAACEtbg4acgQ3ynKj5EtAAAQEaZNk+65x3eKsqtQ2TLG1DLGfGqM+SHvc83D3K+BMeYTY8wKY8xyY0yjijwuAACIPfPnS1OmSLt2+U5SNhUd2bpb0mfW2pMlfZZ3XJTXJY221jaX1EHSpgo+LgAAiDF/+Yv03XdS1aq+k5RNRcvW+ZJey7v8mqQLDr6DMaaFpARr7aeSZK3NsNburuDjAgCAGHPEEVJ8vJSTI+2OoCZR0bJVx1q7Me/yr5LqFHGfppK2G2P+Y4xZZIwZbYyJr+DjAgCAGLR7t9S8ufToo76TlF6J70Y0xkyXdFwRN43Y/8Baa40x9jCP0UlSW0lrJU2SNETSy0U81nWSrpOkBg0alBQNAADEmCpVpEsukU47zXeS0iuxbFlrexzuNmPMb8aYutbajcaYuip6LdY6SYuttT/lfc27kk5XEWXLWvuipBclKTU1tajiBgAAYtwjj/hOUDYVnUacKmlw3uXBkt4r4j7zJdUwxtTOO+4maXkFHxcAAMSw3bulCRPc+q1wV9Gy9ZiknsaYHyT1yDuWMSbVGPNPSbLW5ki6U9JnxpilkoykKDitJAAA8GXaNGnoUGnGDN9JSmasDc/ZutTUVJuWluY7BgAACEPZ2dJXX0kdO0rG+E4jGWMWWGtTi7qN0/UAAICIk5DgTk4dCThdDwAAiFhPPSXdfLPvFMWjbAEAgIi1aZO0fn14L5RnGhEAAESsRx+V4sJ86CjM4wEAABxeftH6/XcpI8NvlsOhbAEAgIi2bp10wgnSiy/6TlI0yhYAAIho9etLo0ZJ557rO0nRWLMFAAAi3vDhvhMcHiNbAAAgKqxaJT3+uBRu+7VTtgAAQFT45BNp5Ejpxx99JzkQZQsAAESFIUOkNWukJk18JzkQa7YAAEBUqFLFfUhuKjEczpkoMbIFAACiSG6udPnl0l/+4jtJIcoWAACIGnFx0jHHSDVq+E5SiGlEAAAQVZ591neCAzGyBQAAolJaWnicoJqyBQAAos706VL79tJ77/lOQtkCAABR6OyzpRdekHr29J2ENVsAACAKxcdL113nO4XDyBYAAIhaH30kjR/vNwMjWwAAIGode6zUvbvfDJQtAAAQtU491XcCphEBAACCirIFAAAQRJQtAACAIKJsAQAABBFlCwAAIIgoWwAAAEFE2QIAAAgiyhYAAEAQUbYAAACCiLIFAAAQRJQtAACAIKJsAQAABBFlCwAAIIgoWwAAAEFE2QIAAAgiyhYAAEAQUbYAAACCiLIFAAAQRJQtAACAIKJsAQAABBFlCwAAIIgoWwAAAEFE2QIAAAgiY631naFIxpjNktb4zhFhjpH0u+8QOACvSfjhNQlPvC7hh9ekbBpaa2sXdUPYli2UnTEmzVqb6jsHCvGahB9ek/DE6xJ+eE0Ch2lEAACAIKJsAQAABBFlK7q86DsADsFrEn54TcITr0v44TUJENZsAQAABBEjWwAAAEFE2YowxphzjDErjTGrjDF3H+Y+lxpjlhtjvjXG/CvUGWNRSa+LMaaBMWamMWaRMeYbY0wfHzljiTHmFWPMJmPMssPcbowxz+S9Zt8YY04NdcZYU4rXZGDea7HUGPOFMSYl1BljUUmvy373a2+MyTbGXByqbNGCshVBjDHxkp6TdK6kFpKuMMa0OOg+J0u6R9KZ1tqWkm4Ldc5YU5rXRdJ9kt621raVdLmkf4Q2ZUyaIOmcYm4/V9LJeR/XSRofgkyxboKKf01+ltTFWtta0kNizVCoTFDxr0v+77nHJX0SikDRhrIVWTpIWmWt/clau0/SW5LOP+g+10p6zlq7TZKstZtCnDEWleZ1sZKOyrtcXdKGEOaLSdbazyVtLeYu50t63TpfSqphjKkbmnSxqaTXxFr7Rf7vLklfSqofkmAxrhT/rUjSLZL+LYn/p5QDZSuy1JP0y37H6/Ku219TSU2NMf8zxnxpjCn2rxUERGlel1GSrjLGrJM0Te4XF/wqzesGf66R9F/fISAZY+pJGiBGf8uNshV9EuSmRbpKukLSS8aYGj4DQZJ7LSZYa+tL6iPpDWMM//0BRTDGnC1Xtob7zgJJ0lhJw621ub6DRKoE3wFQJuslnbDfcf286/a3TtJX1tosST8bY76XK1/zQxMxJpXmdblGeWsirLXzjDFJcucdY0jen9K8bggxY0yypH9KOtdau8V3HkiSUiW9ZYyR3O+tPsaYbGvtu15TRRD+so4s8yWdbIxpbIypJLfQeupB93lXblRLxphj5KYVfwphxlhUmtdlraTukmSMaS4pSdLmkKbEwaZKGpT3rsTTJaVbazf6DhXLjDENJP1H0h+std/7zgPHWtvYWtvIWttI0hRJN1K0yoaRrQhirc02xtws6WNJ8ZJesdZ+a4x5UFKatXZq3m29jDHLJeVIuou/DoOrlK/LHXJTusPkFssPsewoHFTGmP+T+8PjmLy1cg9ISpQka+3zcmvn+khaJWm3pKF+ksaOUrwmIyUdLekfeaMo2ZwIOfhK8bqggthBHgAAIIiYRgQAAAgiyhYAAEAQUbYAAACCiLIFAAAQRJQtAACAIKJsAQAABBFlCwAAIIgoWwAAAEH0/94x4BvDIRQxAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(.5, 1.5, 100)\n",
    "# y = 1 * coef[\"x0\"] + x * coef[\"x\"] + x**2 * coef[\"x2\"]# + x**3 * coef[\"x3\"]\n",
    "\n",
    "y = 1 * coef[0] + np.cos(x) * coef[1] + np.sin(x) * coef[2]# + np.cos(2*x) * coef[3] + np.sin(2*x) * coef[4]\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.plot(x,y, linestyle=\":\", color=\"blue\")\n",
    "\n",
    "plt.scatter(X, Y, color=\"red\", label=\"Y (discounted exercise later)\")\n",
    "\n",
    "plt.scatter(X, continuation, label=\"continuation\", marker=\"x\", color=\"blue\")\n",
    "plt.scatter(X, exercise, label=\"exercise now\", marker=\"+\", color=\"green\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So on this plot we now have the full picture.\n",
    "\n",
    "The dashed line is the estimated function. If we use more terms in the basis function, it gets a better fit. This rough fit is good for illustration purpose\n",
    "\n",
    "Again, the red dots are what the final option value was given the current stock price.\n",
    "\n",
    "The green dots are what we get if we exercise the option now.\n",
    "\n",
    "The blue dots are the theoretical continuation value. It is computed from the fitted continuation taking all values into account. So it is the **expected continuation value**.\n",
    "\n",
    "We can see in some cases, (e.g. the first path (i.e. the first 3 dots from the left)), the exercise now value is very high, it is higher both than the actual value we got later (red dot) and the estimated value (blue dot). It is thus definitely worth it to exercise it.\n",
    "\n",
    "For the last path, so the last 3 dots, we have a different picture. The blue (expected value) is highest, below that is the green (exercise now) and even below the red (what really happened). Anyway, just looking at the green and blue one, it thus makes more sense to exercise now.\n",
    "\n",
    "We can see that this picture, i.e. this data, allows us to determine the optimal choice for each path: is it better to exercise now, or to continue\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### timestamp 1\n",
    "\n",
    "We now do the same thing again. Y is the realised cashflow. So for the cases where we decided to continue (i.e. continuation > exercise), Y will be 0 (since no cash realised). In the cases where we decided to exercise, Y will be the exercise profit.\n",
    "\n",
    "X is again the stock price at current time, which is time t=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1     True\n",
      "2     True\n",
      "3     True\n",
      "4    False\n",
      "5     True\n",
      "6    False\n",
      "7    False\n",
      "8     True\n",
      "dtype: bool\n"
     ]
    }
   ],
   "source": [
    "continued = continuation > exercise\n",
    "\n",
    "# we also consider those that were OTM to be continued, since early exercise would have made no sense\n",
    "# hence we can reindex and fill with True\n",
    "\n",
    "continued = continued.reindex(df.index).fillna(True)\n",
    "\n",
    "print(continued)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    0.000000\n",
       "2    0.000000\n",
       "3    0.000000\n",
       "4    0.122429\n",
       "5    0.000000\n",
       "6    0.310782\n",
       "7    0.244859\n",
       "8    0.000000\n",
       "Name: 2, dtype: float64"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = (k - df[N-1]).map(lambda v: max(v, 0))\n",
    "\n",
    "Y[continued] = 0  # we take realised cash flows, and if we continue no cash is realised\n",
    "\n",
    "discount = np.exp(-mu * 1) # discount again 1 since we realised those next iteration\n",
    "\n",
    "Y = Y * discount\n",
    "Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "X is again computed same as before, X is the stock value at the current moment, only for those that are currently ITM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    1.09\n",
       "4    0.93\n",
       "6    0.76\n",
       "7    0.92\n",
       "8    0.88\n",
       "Name: 1, dtype: float64"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = df[N-2]\n",
    "\n",
    "ITM = X < k\n",
    "\n",
    "X = X[ITM]\n",
    "\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    0.000000\n",
       "4    0.122429\n",
       "6    0.310782\n",
       "7    0.244859\n",
       "8    0.000000\n",
       "Name: 2, dtype: float64"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = Y[ITM]\n",
    "Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have X and Y defined again and can do another regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd7b27a2d00>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAI/CAYAAABAoBw9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABKb0lEQVR4nO3dd5xU1f3/8ddhQRBBEcUSkaKC0kEWLBFBEBsGS9SoREUTTbFFv1+jRkVi9BeNxhq/8atGjQYrRsUSJRYQS5QVKyCKiIj6VYr0zp7fH5eFBRbYZWf2zuy8no/HPmbOnbtzP+xQ3pxz7jkhxogkSZKyo07aBUiSJNVmhi1JkqQsMmxJkiRlkWFLkiQpiwxbkiRJWWTYkiRJyqK6aRewIdtvv31s1apV2mVIkiRt0jvvvDMzxtisotdyNmy1atWKkpKStMuQJEnapBDCFxt6zWFESZKkLDJsSZIkZZFhS5IkKYtyds6WJKn2Wb58OdOnT2fJkiVplyJtlgYNGtC8eXPq1atX6e8xbEmSasz06dNp3LgxrVq1IoSQdjlSlcQYmTVrFtOnT6d169aV/j6HESVJNWbJkiVst912Bi3lpRAC2223XZV7Zg1bkqQaZdBSPtuc37+GLUmSpCwybEmSCkaMkQMOOIB//etfq4899thjHHbYYeud17dvX+bNm7feewwdOpQbbrgBgCFDhvDiiy9mt+iNeO+993juueeq/H19+vRJfeHwTPzsRo0axZFHHrnRczb3ZzRjxoz1fl9sLsOWJKlghBC44447uPDCC1myZAkLFizgd7/7Hbfffvta5z333HN06dKFrbfeeqPvd9VVV3HwwQdns+SN2twgUVNijJSWllb4Wk397DbnZ7RixQqaNWvGzjvvzOuvv17tGgxbkqTcNWwYtGoFdeokj8OGVfstO3bsyI9+9COuu+46rrrqKk499VR23333dS47jKOOOmp1+5prrqFt27YccMABTJo0afXxwYMHM3z4cAAuueQS2rdvT+fOnfnv//5vAL799luOOeYYunTpQpcuXXjjjTcAuPHGG+nYsSMdO3bk5ptvBmDq1Kl07Nhx9XvfcMMNDB06FEh6oi6++GJ69uxJ27ZtGTNmDMuWLWPIkCE88sgjdO3alUceeYSFCxdyxhln0LNnT7p168ZTTz0FwOLFiznxxBNp164dxxxzDIsXL67wZ/POO+/Qu3dvunfvzqGHHso333zD3Llz2XPPPVf/uk866STuuusuAK6//np69OhB586dufLKK1f/Ovbcc09OPfVUOnbsyJdffsl1111Hp06d6NKlC5dcckmlfnYzZszgxz/+MT169KBHjx6bDD1vv/02++23H926dWP//fdn0qRJVfoZ3XfffQwcOJC+ffvSr18/AI4++miGZeD3HDHGnPzq3r17lCTVLhMmTKj8yf/4R4wNG8YIa74aNkyOV9OCBQti27ZtY8eOHeOSJUvWe71FixZx3rx5McYYS0pKYseOHePChQvj3Llz4+677x6vv/76GGOMp512WnzsscfizJkzY9u2bWNpaWmMMcbvv/8+xhjjCSecEG+66aYYY4wrVqyIc+bMWf1+CxYsiPPnz4/t27eP48aNi59//nns0KHD6hquv/76eOWVV8YYY+zdu3e88MILY4wxPvvss7Ffv34xxhjvvffeePbZZ6/+nksvvTQ+8MADq2to06ZNXLBgQfzzn/8cTz/99BhjjO+//34sKiqKY8eOXevXvGzZsrjffvvF7777LsYY48MPP7z6e0aOHBn33Xff+NBDD8VDDz00xhjjCy+8EM8888xYWloaV65cGQcMGBBHjx4dP//88xhCiG+++WaMMcbnnnsu7rfffnHhwoUxxhhnzZpVqZ/dSSedFMeMGRNjjPGLL76Ie+2113qf0yuvvBIHDBgQY4xx7ty5cfny5THGGP/973/HY489tko/o3vvvTfusssuq+uLMcbp06fHjh07rnfdin4fAyVxA5nGdbYkSbnpsstg0aK1jy1alBwfNKhab73VVlvxk5/8hEaNGlG/fv31Xp89ezaNGzcGYMyYMRxzzDE0bNgQgIEDB653/jbbbEODBg342c9+xpFHHrl6HtHLL7/M/fffD0BRURHbbLMNr732GscccwxbbbUVAMceeyxjxoyp8H3LO/bYYwHo3r07U6dOrfCckSNHMmLEiNVzypYsWcK0adN49dVXOe+88wDo3LkznTt3Xu97J02axEcffUT//v0BWLlyJTvvvDMA/fv357HHHuPss8/m/fffX32tkSNH0q1bNwAWLFjAp59+SosWLWjZsiX77rsvAC+++CKnn3766p9f06ZNK/Wze/HFF5kwYcLq8+bNm8eCBQto1KhRhb/2uXPnctppp/Hpp58SQmD58uVV+hmV/TrL17fDDjvw9ddfV/g+VWHYkiTlplX/AFb6eBXVqVOHOnUqnk1Tt25dSktLN/h6Ree//fbbvPTSSwwfPpy//OUvvPzyy1Wqp+yaZdZdy6ksFBYVFbFixYoK3yPGyOOPP86ee+5ZpWuXfW+HDh14880313uttLSUiRMn0rBhQ77//nuaN29OjJFLL72UX/ziF2udO3Xq1NVBsjI29LMrLS3lP//5Dw0aNKjU+1xxxRUcdNBBPPHEE0ydOpU+ffps8NdZ0c/orbfeWq/uJUuWsOWWW1b617IhztmSJOWmFi2qdjyD9txzT6ZMmQLAgQceyJNPPsnixYuZP38+Tz/99HrnL1iwgLlz53LEEUdw0003re796devH3/961+BpKdo7ty59OrViyeffJJFixaxcOFCnnjiCXr16sWOO+7Id999x6xZs1i6dCnPPPPMJuts3Lgx8+fPX90+9NBDue2220hGteDdd99d/Wt48MEHAfjoo4/44IMPKvw1z5gxY3XYWr58OePHjwfgpptuol27djz44IOcfvrpLF++nEMPPZR77rmHBQsWAPDVV1/x3Xffrfe+/fv3595772XRql7K2bNnV+pnd8ghh3DbbbetPu+9997b6M9i7ty57LLLLkAy/6qqP6OKfPLJJ2vNo9tchi1JUm665hpYNfS0WsOGyfEsGzBgAKNGjQJg77335ic/+QldunTh8MMPp0ePHuudP3/+fI488kg6d+7MAQccwI033gjALbfcwiuvvEKnTp3o3r07EyZMYO+992bw4MH07NmTffbZh5///Od069aNevXqMWTIEHr27En//v3Za6+9NlnnQQcdxIQJE1ZP/r7iiitYvnw5nTt3pkOHDlxxxRUA/OpXv2LBggW0a9eOIUOG0L179/Xea4sttmD48OFcfPHFdOnSha5du/LGG28wadIk7r77bv785z/Tq1cvDjzwQK6++moOOeQQTj75ZPbbbz86derEcccdt1aoKXPYYYcxcOBAiouL6dq16+rhu0397G699VZKSkro3Lkz7du354477tjoz+K3v/0tl156Kd26dVur56+yP6OKvPLKKwwYMGCj162MUJbsck1xcXFMew0QSVJmTZw4kXbt2lX+G4YNS+ZoTZuW9Ghdc02152tVxjfffMOpp57Kv//976xfS7nrwAMP5KmnnmLbbbdd63hFv49DCO/EGIsrep+MzNkKIRwG3AIUAXfHGK9d5/WbgINWNRsCO8QYm2Ti2pKkWmzQoBoJV+vaeeedOfPMM5k3b94m19pS7TRjxgwuvPDC9YLW5qh22AohFAG3A/2B6cDYEMKIGOPqWwhijBeUO/9coFt1rytJUjadcMIJaZegFDVr1oyjjz46I++ViTlbPYHJMcYpMcZlwMPAURs5/yTgoQxcV5IkKedlImztAnxZrj191bH1hBBaAq2Bqt0PK0mSlKdq+m7EE4HhMcaVFb0YQjgrhFASQiiZMWNGDZcmSZKUeZkIW18Bu5ZrN191rCInspEhxBjjnTHG4hhjcbNmzTJQ2sYtWgQVLDUiSdIGTZ06dfWaVQAlJSWrV2fPpCeffHKtFdSHDBnCiy++mPHrKPsyEbbGAm1CCK1DCFuQBKoR654UQtgL2BZYf2nalAwaBAMGwAZW9JckaT3rhq3i4mJuvfXWjF9n3bB11VVXcfDBB2f8Osq+aoetGOMK4BzgBWAi8GiMcXwI4aoQQvmNnk4EHo45tLDXJZckS7jUddMiScpJ6/6LkYl/Qe6//346d+5Mly5dOOWUU5g6dSp9+/alc+fO9OvXb/U+eYMHD+a8885j//33Z7fddmP48OEAXHLJJYwZM4auXbty0003MWrUqNX7+Q0dOpQzzjiDPn36sNtuu60OYVOnTl1rJfIbbriBoUOHAnDXXXfRo0cPunTpwo9//GMWLVrEG2+8wYgRI7jooovo2rUrn332GYMHD15dw0svvUS3bt3o1KkTZ5xxBkuXLgWgVatWXHnlley999506tSJjz/+uPo/MFVbRuZsxRifizG2jTHuHmO8ZtWxITHGEeXOGRpjvCQT18uUffaBAw+EENKuRJK0rqFD4YIL1gSsGJP2qoyyWcaPH8/VV1/Nyy+/zPvvv88tt9zCueeey2mnncYHH3zAoEGD1hoS/Oabb3jttdd45plnuOSS5J+wa6+9ll69evHee+9xwQUXrHeNjz/+mBdeeIG3336b3//+9xvcELnMsccey9ixY3n//fdp164df/vb39h///0ZOHAg119/Pe+99x6777776vOXLFnC4MGDeeSRR/jwww9ZsWLF6i2BALbffnvGjRvHr371q/VWa1c6Cn67ngUL4Kqr4PXX065EklQmRpgzB265ZU3guuCCpD1nzub3cL388sscf/zxbL/99gA0bdqUN998k5NPPhmAU045hddee231+UcffTR16tShffv2fPvtt5W6xoABA6hfvz7bb789O+ywwya/76OPPqJXr1506tSJYcOGrd6PcEMmTZpE69atadu2LQCnnXYar7766urXjz32WAC6d+/O1KlTK1WzsqvgB9CKiuAvf0l6t374w7SrkSRB8nfyTTclz2+5JfkCOP/85HhNjUjUr19/9fPKzoIp/z1FRUWsWLGCunXrUlpauvr4kiVLVj8fPHgwTz75JF26dOG+++5bvSdjdWsuu7bSV/A9W1tuCZ98AhvZh1KSlILygatMdYNW3759eeyxx5g1axYAs2fPZv/99+fhhx8GYNiwYfTq1Wuj79G4ceMKN1zemB133JHvvvuOWbNmsXTpUp555pnVr82fP5+dd96Z5cuXM2zYsE1eZ88992Tq1KlMnjwZgAceeIDevXtXqR7VrIIPWwBNmiSP5f6jIUlKWdnQYXnl53Btjg4dOnDZZZfRu3dvunTpwoUXXshtt93GvffeS+fOnXnggQe4pawbbQM6d+5MUVERXbp04aZ10+AG1KtXjyFDhtCzZ0/69+/PXnvttfq1P/zhD+yzzz788Ic/XOv4iSeeyPXXX0+3bt347LPPVh9v0KAB9957L8cffzydOnWiTp06/PKXv6ziT0I1KeTQzYFrKS4ujiUlJTV2vYcegnPOgY8/hhpY4kuSCtLEiRNp167dJs8rP0erbOhw3bY3NyktFf0+DiG8E2Msruj8gp+zVaZrVzj6aFi2LO1KJEkhJKMO5YNVWSdSkyYGLeUXw9Yq7drB3/6WdhWSpDJDhyY9XGXBqixwGbSUb5yztY7PPnMZCEnKFesGK4OW8pE9W+s48cRk+5733ku7EkmSVBsYttZx552w005pVyFJkmoLw9Y6unVLuwJJklSbOGerAl99BccfD2PHpl2JJCmfHXHEEcyZMyftMpQyw1YFtt4aSkqSyfKSJG3MxrbEee6552hStnK2CpZhqwKNG8PkyclkeUlS7fKPf/yDnj170rVrV37xi1+wcuVKxo4dS+fOnVmyZAkLFy6kQ4cOfPTRRyxcuJAzzjiDnj170q1bN5566ikA7rvvPgYOHEjfvn3p168fCxYs4PTTT6dTp0507tyZxx9/HIBWrVoxc+ZMFi5cyIABA+jSpQsdO3bkkUceAeCdd96hd+/edO/enUMPPZRvvvlmvXoHDx7Meeedx/77789uu+3G8OHDgWSvxosuuoiOHTvSqVOn1e959tlnM2LECACOOeYYzjjjDADuueceLrvssuz+cFUh52xtQFFR8jh9OjRvnm4tklTI+tzXB4BRg0dV+70mTpzII488wuuvv069evX49a9/zbBhwzj11FMZOHAgl19+OYsXL+anP/0pHTt25He/+x19+/blnnvuYc6cOfTs2ZODDz4YgHHjxvHBBx/QtGlTLr74YrbZZhs+/PBDAL7//vu1rvv888/zgx/8gGeffRaAuXPnsnz5cs4991yeeuopmjVrxiOPPMJll13GPffcs17d33zzDa+99hoff/wxAwcO5LjjjuOf//wn7733Hu+//z4zZ86kR48eHHjggfTq1YsxY8YwcOBAvvrqq9UBbsyYMZxoL0IqDFsb8cc/wtVXw5dfQtOmaVcjSaqul156iXfeeYcePXoAsHjxYnbYYQcAhgwZQo8ePWjQoAG33norACNHjmTEiBHccMMNACxZsoRp06YB0L9/f5qu+sfhxRdfXL2ZNcC222671nU7derEf/3Xf3HxxRdz5JFH0qtXLz766CM++ugj+vfvD8DKlSvZeeedK6z76KOPpk6dOrRv355vv/0WgNdee42TTjqJoqIidtxxR3r37s3YsWPp1asXN998MxMmTKB9+/Z8//33fPPNN7z55purf12qWYatjTjqKGjQAOrXT7sSSSo8ZT1ao78YvVa7Oj1cMUZOO+00/vjHP6732qxZs1iwYAHLly9nyZIlbLXVVsQYefzxx9lzzz3XOvett95iq622qvR127Zty7hx43juuee4/PLL6devH8cccwwdOnTgzTff3OT31y/3D9Gm9jTeZZddmDNnDs8//zwHHnggs2fP5tFHH6VRo0Y0bty40jUrc5yztRHt2ycbn1bhz5MkKYf169eP4cOH89133wEwe/ZsvvjiCwB+8Ytf8Ic//IFBgwZx8cUXA3DooYdy2223rQ447777boXv279/f26//fbV7XWHEb/++msaNmzIT3/6Uy666CLGjRvHnnvuyYwZM1aHreXLlzN+/PhK/1p69erFI488wsqVK5kxYwavvvoqPXv2BGDffffl5ptvXj2seMMNN9CrV69Kv7cyy56tTYgRnn0W6tSBI45IuxpJKhxlPViZnLPVvn17rr76ag455BBKS0upV68et99+O6NHj6ZevXqcfPLJrFy5kv3335+XX36ZK664gt/85jd07tyZ0tJSWrduzTPPPLPe+15++eWcffbZdOzYkaKiIq688kqOPfbY1a9/+OGHXHTRRdSpU4d69erx17/+lS222ILhw4dz3nnnMXfuXFasWMFvfvMbOnToUKlfyzHHHMObb75Jly5dCCHwpz/9iZ1Wrcrdq1cvRo4cyR577EHLli2ZPXu2YStFYVPdkWkpLi6OJSUlaZdBjNC9O2y/PYwcmXY1kpTfJk6cSLt27ar0PZkMW1ImVPT7OITwToyxuKLz7dnahBDgiSfgBz9IuxJJKkyGLOU7w1YltGyZPJaWJsOJkiRJlWV0qKQJE6BDB3jjjbQrkSRJ+cSwVUktW8LOO8Py5WlXIkn5LVfnCkuVsTm/fx1GrKSttoKXX067CknKbw0aNGDWrFlst912hBDSLkeqkhgjs2bNokGDBlX6PsNWFS1dCm+/Dd5BK0lV17x5c6ZPn86MGTPSLkXaLA0aNKB5FffxM2xV0eWXw623wrRpsOOOaVcjSfmlXr16tG7dOu0ypBpl2Kqis8+GQw+FVVtpSZIkbZRhq4patUq+JEmSKsO7ETdDaSn88Y/wt7+lXYkkScp1hq3NUKdOsnXP66+nXYkkScp1DiNupueegy23TLsKSZKU6+zZ2kxlQev775NhRUmSpIoYtqph3DjYdVd45pm0K5EkSbnKsFUNnTvDGWdA27ZpVyJJknKVc7aqoW7dZIFTSZKkDbFnKwOmTYM770y7CkmSlIsMWxlw331wzjnwzTdpVyJJknKNYSsDzjsPJk+GnXdOuxJJkpRrnLOVAU2aJF8AMUIIaVYjSZJyiT1bGRIjDB4MF12UdiWSJCmXGLYyJARo1Ai22irtSiRJUi5xGDGD/vKXtCuQJEm5xp6tLBg3DhYsSLsKSZKUCwxbGTZhAnTvDnffnXYlkiQpFziMmGHt28P998PAgWlXIkmScoFhKwtOOSXtCiRJUq5wGDFL3nwTTj8dVq5MuxJJkpQmw1aWTJ8OI0fC55+nXYkkSUqTYStLjj0WpkyBPfZIuxJJkpQmw1aWFBVB/frJyvKzZqVdjSRJSktGwlYI4bAQwqQQwuQQwiUbOOeEEMKEEML4EMKDmbhuPjjySDj++LSrkCRJaan23YghhCLgdqA/MB0YG0IYEWOcUO6cNsClwA9jjN+HEHao7nXzxcknw4oVblAtSVKhysTSDz2ByTHGKQAhhIeBo4AJ5c45E7g9xvg9QIzxuwxcNy8MGpR2BZIkKU2ZGEbcBfiyXHv6qmPltQXahhBeDyH8J4RwWAaumzeWL4cHHkhWl5ckSYWlpibI1wXaAH2Ak4C7QghN1j0phHBWCKEkhFAyY8aMGiot++bPh1//Gu67L+1KJElSTcvEMOJXwK7l2s1XHStvOvBWjHE58HkI4ROS8DW2/EkxxjuBOwGKi4tjBmrLCU2bwtixsOeeaVciSZJqWiZ6tsYCbUIIrUMIWwAnAiPWOedJkl4tQgjbkwwrTsnAtfPGXnslE+RLS9OuRJIk1aRqh60Y4wrgHOAFYCLwaIxxfAjhqhBC2XbMLwCzQggTgFeAi2KMBbf61MsvQ5s2yerykiSpMGRkI+oY43PAc+scG1LueQQuXPVVsHbbDVq3hnnz0q5EkiTVlIyELVVOq1bw4otpVyFJkmqS2/WkYN48ePXVtKuQJEk1wbCVgrPPhoEDYdGitCuRJEnZZthKwWWXwUsvQcOGaVciSZKyzTlbKdhrr7QrkCRJNcWerZSsWAG//S387/+mXYkkScome7ZSUrcujBuXhC5JklR7GbZS9PzzSeiSJEm1l8OIKSoLWl9+CcuXp1uLJEnKDsNWyj78EHbfHYYNS7sSSZKUDYatlHXsCFdeCX37pl2JJEnKBmcMpSyEZN0tSZJUO9mzlSM+/hguvhhKS9OuRJIkZZJhK0eUlMDttyehS5Ik1R6GrRxx4okwdSq0b592JZIkKZMMWzmibl3YfvvkuRtUS5JUexi2cswvfwkHHwwxpl2JJEnKBO9GzDG9ekHLlrBypavLS5JUG/jPeY4ZNCjtCiRJUiY5jJiDYoSXXoJ33km7EkmSVF2GrRy0bBmccgr86U9pVyJJkqrLYcQcVL8+PP88tG2bdiWSJKm6DFs5qnPn5DHGZEsfSZKUnxxGzGEffwzdusHbb6ddiSRJ2lyGrRzWvDk0bAgLFqRdiSRJ2lwOI+awRo3gjTfSrkKSJFWHPVt5YMUK+Pe/065CkiRtDsNWHvjLX+CQQ+CDD9KuRJIkVZXDiHngjDNg992hU6e0K5EkSVVl2MoDW28NP/pR2lVIkqTN4TBiHvnb3+BXv0q7CkmSVBWGrTwyfXqy9tbSpWlXIkmSKsuwlUcuuwxeeSXZzkeSJOUHw1Yeqbtqht3cufDVV+nWIkmSKsewlWdWroSuXeG889KuRJIkVYZ3I+aZoiL44x+hTZu0K5EkSZVh2MpDJ56YdgWSJKmyHEbMU3PmwEUXwbvvpl2JJEnaGMNWngoB7r0XxoxJuxJJkrQxDiPmqW22gSlTktXlJUlS7rJnK4+VBa3vvku3DkmStGGGrTz33HPQvDmMHZt2JZIkqSKGrTzXqxeccw7sskvalUiSpIo4ZyvPNW4MN96YdhWSJGlD7NmqJSZOhGuvTbsKSZK0LsNWLTFiBFxzjXsmSpKUawxbtcS558Lnnzt3S5KkXGPYqiUaNoTtt0+eL1qUbi2SJGkNw1Yt8/Ofw2GHQYxpVyJJksC7EWudAw6APfaA0lIoKkq7GkmSZNiqZQYPTrsCSZJUnsOItVCM8OKL8PrraVciSZIyErZCCIeFECaFECaHEC6p4PXBIYQZIYT3Vn39PBPXVcVWroSzzoLrrku7EkmSVO1hxBBCEXA70B+YDowNIYyIMU5Y59RHYoznVPd62rS6deHZZ2G33dKuRJIkZaJnqycwOcY4Jca4DHgYOCoD76tqaNcO6tdPJsqXlqZdjSRJhSsTYWsX4Mty7emrjq3rxyGED0IIw0MIu2bgutqEr76Cbt3g0UfTrkSSpMJVUxPknwZaxRg7A/8G/l7RSSGEs0IIJSGEkhkzZtRQabXXzjsnQ4mNG6ddiSRJhSsTSz98BZTvqWq+6thqMcZZ5Zp3A3+q6I1ijHcCdwIUFxe7LGc11akDTzyRdhWSJBW2TPRsjQXahBBahxC2AE4ERpQ/IYSwc7nmQGBiBq6rSlqxAv7xD1i6NO1KJEkqPNXu2YoxrgghnAO8ABQB98QYx4cQrgJKYowjgPNCCAOBFcBsYHB1r6vKe/VVOOWUpKfr5JPTrkaSpMISYo5uoldcXBxLSkrSLqNWiBFGj4bevSGEtKuRJKn2CSG8E2Msrug1V5AvACFAnz7JY45ma0mSai3DVgF5+mno3BnmzUu7EkmSCodhq4DstBNstx3MmrXpcyVJUmZkYukH5YkePWDUqLSrkCSpsNizVYDmzoWRI9OuQpKkwmDYKkC//S0cfTTMmZN2JZIk1X6GrQL0u9/Ba69BkyZpVyJJUu3nnK0C1LJl8iVJkrLPnq0CNnQonHlm2lVIklS7GbYK2PLlyVdpadqVSJJUezmMWMCuvtrteyRJyjZ7tgpYWdD6/HOYODHdWiRJqq0MWwVu5cpkg+oLLki7EkmSaieHEQtcURH8/e/Qtm3alUiSVDsZtsRBB6VdgSRJtZfDiAKS1eRPOAEeeSTtSiRJql0MWwKgcWP46iuYNSvtSiRJql0cRhSQzN167TWXgpAkKdPs2dJqZUHrzTdhwYJ0a5EkqbYwbGktH38M++8Pf/1r2pVIklQ7OIyotey1VzJJfsCAtCuRJKl2MGxpPSeckHYFkiTVHg4jqkIffpisLD9tWtqVSJKU3wxbqtA228DXX8MXX6RdiSRJ+c1hRFWoRQuYNAnqGMclSaoW/ynVBtWpA6Wl8NJLaVciSVL+Mmxpo+69Fw4+GN54I+1KJEnKTw4jaqMGDYJGjWDffdOuRJKk/GTY0kY1aAA/+UnaVUiSlL8cRlSlPP009O0Ly5alXYkkSfnFsKVKKSpK9kv89tu0K5EkKb84jKhKOfzw5Ktss2pJklQ59mypUkJIvhYtgtGj065GkqT8YdhSlfzXf8ERR8Ds2WlXIklSfjBsqUouvhheeAGaNk27EkmS8oNztlQlrVolXwAxOodLkqRNsWdLm+XGG+GEE9KuQpKk3GfY0maJMXlcsiTdOiRJynUOI2qzXHihQ4iSJFWGPVvaLGVB6+uv4fXX061FkqRcZs+WquXkk2H6dPjkE6hjdJckaT2GLVXLrbdCo0YGLUmSNsSwpWrp3HnNc5eCkCRpffZHqNpihLPOgnPOSbsSSZJyjz1bqrYQoEkTqFvX3i1JktZl2FJG/OlPaVcgSVJuchhRGTVuHLz1VtpVSJKUO+zZUsasXAk/+Qnsuiu8/HLa1UiSlBsMW8qYoiIYPnzNRtWSJMmwpQzr0iV5jBFWrIB69dKtR5KktDlnSxm3dCn07Qu//33alUiSlL6MhK0QwmEhhEkhhMkhhEs2ct6PQwgxhFCciesqN9WvD506QevWaVciSVL6qj2MGEIoAm4H+gPTgbEhhBExxgnrnNcYOB/wXrUCcOutaVcgSVJuyETPVk9gcoxxSoxxGfAwcFQF5/0BuA5YkoFrKg/ECE88AWPHpl2JJEnpyUTY2gX4slx7+qpjq4UQ9gZ2jTE+m4HrKU8sXgxnn20vlySpsGX9bsQQQh3gRmBwJc49CzgLoEWLFtktTFnXsGGy3tYee6RdiSRJ6clEz9ZXwK7l2s1XHSvTGOgIjAohTAX2BUZUNEk+xnhnjLE4xljcrFmzDJSmtO21V7Jn4rJlsMQBZElSAcpE2BoLtAkhtA4hbAGcCIwoezHGODfGuH2MsVWMsRXwH2BgjLEkA9dWHpg3Dzp2dP9ESVJhqnbYijGuAM4BXgAmAo/GGMeHEK4KIQys7vsr/229Nfz4x7DffmlXIklSzQsxxrRrqFBxcXEsKbHzS5Ik5b4QwjsxxgrXEXUFedWYZcvgxhvhLVdakyQVEPdGVI1ZtgxuuAG++Qb22SftaiRJqhmGLdWYRo1g3DjYaae0K5EkqeY4jKgaVRa0Zs6EBQvSrUWSpJpg2FKNmz0b9twTrr467UokSco+hxFV45o2hSuugEMOSbsSSZKyz7ClVPzmN2lXIElSzXAYUamZPx/OPRdefDHtSiRJyh7DllKzxRYwcmRyh6IkSbWVw4hKTf368P770KBB2pVIkpQ99mwpVWVBa+JE+O67dGuRJCkbDFtK3ezZUFwMQ4emXYkkSZnnMKJS17Qp3H8/HHhg2pVIkpR5hi3lhB//eM3zGCGE9GqRJCmTHEZUzpg3D448Eu69N+1KJEnKHMOWckbjxlBaCitXpl2JJEmZ4zCickYI8OyzDiFKkmoXe7aUU8qC1r//DePHp1uLJEmZYNhSzlmwAE46Ca67Lu1KJEmqPocRlXMaNUq28enQIe1KJEmqPsOWctLeeyePK1bA8uWw5Zbp1iNJ0uZyGFE5a8kS6NEDLr887UokSdp89mwpZzVoAEcdBV27pl2JJEmbz7ClnOZ+iZKkfOcwonJeaSnceSc8/njalUiSVHWGLeW8GOHuu+Gxx9KuRJKkqnMYUTmvqAj+9S9o2jTtSiRJqjp7tpQXttsuWV1+zhz4/PO0q5EkqfLs2VLeiBF69042rB4zxj0UJUn5wbClvBECXHst7LCDQUuSlD8MW8orhx++5nmMhi5JUu5zzpby0pVXwq9+lXYVkiRtmmFLeWnZsuRr5cq0K5EkaeMcRlRe+n//zyFESVJ+sGdLeaksaH3+OTzxRLq1SJK0MYYt5bVLLknmbi1ZUv336nNfH/rc16f6byRJUjkOIyqv3XhjsndigwZpVyJJUsUMW8pru+yy5vmsWclK81VV1ps1+ovRa7VHDR5VveIkScJhRNUSQ4ZAly4wf37alUiStDZ7tlQrDBiQbFi9xRZV/96yHix7tCRJ2WDYUq2wzz7JlyRJucawpVrl7bfhL3+Be+9Nerqqwh4tSVI2OGdLtcrUqfDKK8mjJEm5wLClWuX44+Hjj2H33dOuRJKkhGFLtUoIsNVWydpbzzwDMaZdkSSp0Bm2VCs99BD86Efw0ktpVyJJKnSGLdVKP/kJPP449OtXyW8YNgxatYI6dZLHYcOyWJ0kqZB4N6Jqpbp14dhjk+dLlmxiO59hw+Css2DRoqT9xRdJG2DQoKzWKUmq/ezZUq323nuw227w6qsbOemyy9YErTKLFiXHJUmqJsOWarW2bWG//WCbbTZy0rRpVTsuSVIVOIyoWq1hw2Tu1ka1aJEMHVZ0XJKkarJnSwVh8WL43e/ggw8qePGaa5JUVl7DhslxSZKqKSNhK4RwWAhhUghhcgjhkgpe/2UI4cMQwnshhNdCCO0zcV2pshYtgr/9Df71rwpeHDQI7rwTWrZMFupq2TJpOzlekpQBIVZz1ccQQhHwCdAfmA6MBU6KMU4od87WMcZ5q54PBH4dYzxsY+9bXFwcS0pKqlWbVN6sWbDddmlXIUmqjUII78QYiyt6LRM9Wz2ByTHGKTHGZcDDwFHlTygLWqtsBbiut2pcWdCaMqXiKVqSJGVDJibI7wJ8Wa49Hdhn3ZNCCGcDFwJbAH0zcF2pypYuhQMOgJ494ckn065GklQIamyCfIzx9hjj7sDFwOUVnRNCOCuEUBJCKJkxY0ZNlaYCUr8+3HMP3HZb2pVIkgpFJsLWV8Cu5drNVx3bkIeBoyt6IcZ4Z4yxOMZY3KxZswyUJq3vsMNg11W/Y5ctS7cWSVLtl4mwNRZoE0JoHULYAjgRGFH+hBBCm3LNAcCnGbiuVC0//3myh2I17xGRJGmjqj1nK8a4IoRwDvACUATcE2McH0K4CiiJMY4AzgkhHAwsB74HTqvudaXq6tgR5s1LwlYIaVcjSaqtqr30Q7a49INq0rqBywAmSaqKbC/9IOW1n/8cevRYM5wYI1xwAQwdmmpZkqRawrClghYjTJgA77wDv/jFmqB1yy0wZ47zuSRJ1edG1CpoIcDo0XDeeXDHHXDXXcnx88+Hm25yKFGSVH3O2ZJIerDqlOvnLS01aEmSKs85W9JGlA0dlnfBBQ4hSpIyw2FEFbTyc7TOOQd69YLXX0/a4FCiJKn6DFsqaCFAkyZrz9E6/nhYvjw5btCSJFWXYUsFb+jQtdfVmjgR/vlP+OtfUy1LklRLOGdLYu0erDZt4NBDYffd06tHklR72LMlraNePfj739OuQpJUW9izJW3A0qVw8cXw2GNpVyJJymeGLWkDiopg1KhkdXlJkjaXw4jSBtStm4StLbdMuxJJUj6zZ0vaiLKg9dln8Pzz6dYiScpP9mxJlXDOOcmSEJ9+mkyglySpsgxbUiXccUeyd6JBS5JUVQ4jSpXQsiXsumvyfMqUdGvJZX3u60Of+/qkXYYk5RTDllQFN98MHTrA5MlpVyJJyhcOI0pV8JOfwKJFSU+X1ijrzRr9xei12qMGj0qnIEnKIYYtqQp23hl+97vkeWlpMo9LkqSNMWxJm+Hjj+H44+Gee6BHj7SrSV9ZD5Y9WpK0Pv9fLm2GHXeExo2TLX0kSdoYe7akzbDttvD66xBC2pXkFnu0JGl99mxJmymEZN7WrbfCs8+mXY0kKVcZtqRqWL4c/vY3ePzxtCuRJOUqw5ZUDfXrw8svJ4Gr4A0bBq1aJbdotmqVtCVJhi2purbbLhlSnDkTRo5Mu5qUDBsGZ50FX3wBMSaPZ51l4JIkDFtSxpx/Ppx4IixYkHYlKbjssmS11/IWLUqOS1KBM2xJGfKnPyVDio0apV1JCqZNq9pxSSoghi0pQ3bZBbp2TZ5PnZpmJSlo0aJqxyWpgBi2pAx7+mnYYw945ZW0K6lB11wDDRuufaxhw+S4JBU4w5aUYf36waWXFtg2PoMGwZ13Jjt0h5A83nlnclySClyIMaZdQ4WKi4tjSUlJ2mVI1bJyZbISgivNS1LtFkJ4J8ZYXNFr9mxJWTJnDvTuDXfckXYlkqQ0GbakLNl6a2jePNlHUZJUuNyIWsqSOnXg4YfTrkKSlDZ7tqQa8MQTcNFFaVchSUqDYUuqAW+/DaNHr7/IuiSp9jNsSTXgqqvgtdfWX4pKklT7GbakGlCvHmyxBSxeDNddB8uWpV2RJKmmGLakGjRqFFxyCfz732lXIkmqKd6NKNWgww+H8eOhffu0K5Ek1RR7tqQaVha0xo+H6dPTrUWSlH2GLSkFixcneyied17alUiSss1hRCkFW24J//gHdOyYdiWSpGyzZ0tKycEHw047QYzwxRdpVyNJyhbDlpSyK6+EvfeG//u/tCuRJGWDw4hSyk45BRo3hh12SLsSSVI22LMlpaxNm2TfxDp1YMmStKuRJGWaYUvKEZ99Bu3awT//mXYlkqRMMmxJOaJ5c+jRA3bdNe1KJEmZlJGwFUI4LIQwKYQwOYRwSQWvXxhCmBBC+CCE8FIIoWUmrivVJvXrw6OPJoELkrsUJUn5r9phK4RQBNwOHA60B04KIay7Gcm7QHGMsTMwHPhTda8r1WZ/+hP8/OcGLkmqDTLRs9UTmBxjnBJjXAY8DBxV/oQY4ysxxkWrmv8BmmfgulKttWgRLFwIy5enXYkkqboysfTDLsCX5drTgX02cv7PgH9l4LpSrTVkCISQfEmS8luNTpAPIfwUKAau38DrZ4UQSkIIJTNmzKjJ0qScUqdOErS+/RZ++lOYOTPtiiRJmysTYesroPz9U81XHVtLCOFg4DJgYIxxaUVvFGO8M8ZYHGMsbtasWQZKk/Lbl1/Cc8/Bu++mXYkkaXNlYhhxLNAmhNCaJGSdCJxc/oQQQjfgf4HDYozfZeCaUkEoLoapU2HrrdOuRJK0uardsxVjXAGcA7wATAQejTGODyFcFUIYuOq064FGwGMhhPdCCCOqe12pUJQFrWefhSefTLUUSdJmCDFH7y0vLi6OJSUlaZch5YTSUjjgAKhXD0aNcuK8JOWaEMI7Mcbiil5zI2opD9SpA088kfRyGbQkKb+4XY+UJ3bcEbbcEpYtgzvuSHq7JEm5z7Al5ZnHH4df/QpeeSXtSiRJleEwopRnTjwRdtsN9tnY0sGSpJxhz5aUZ0JYE7QmToTx49OtR5K0cfZsSXlq5Uo45hjYbjt47TUnzktSrjJsSXmqqAgefBCaNTNoSVIucxhRymN77w277goxJr1bkqTcY9iSaoFHHoFevWDkyLQrkSSty2FEqRY47jhYuBAOPjjtSiRJ67JnS6oF6taFn/0sWWl+9myYNi3tiiRJZezZkmqRGOHww5M7Fd9+OwlfkqR0GbakWiQEuPZaqF/foCVJucKwJdUyBx205vnkybDHHunVIklyzpZUa40eDXvtBY8+mnYlklTYDFtSLfXDH8KQIckcLklSegxbUi1Vt24Stho3huXL4csv065IkgqTYUsqAD/7GfTuDYsWpV2JJBUeJ8hLBeCcc6BPH2jYMO1KJKnwGLakAtCzZ/IFMH06/OAHLg0hSTXFv26lAjJtGnTpAtdckyyAWt66bUlSZhi2pAKy665wwQUwc2byWBawYkzaQ4emWp4k1UqGLamAhACXXZY83nJLMnG+LGjdcgvMmWMPlyRlWog5+jdrcXFxLCkpSbsMqVaKEfr2hVGj1hw7/3y46aYkiEmSqiaE8E6Msbii1+zZkgpQCPDQQ2sfM2hJUnYYtqQCFGOyYXV5v/ylQ4iSlA2GLanAlJ+jdf75sGIF7LAD3Hkn/OY3Bi5JyjTDllRgQoAmTdbM0Soqgv/9XzjmGNh2W4cSJSnTnCAvFagY1w5WZe0PP4SOHQ1dklQVTpCXtJ51w1QI8NFH0L073HxzKiVJUq3kdj2SVuvQAf78ZzjllLQrkaTaw54tSauFAOeem8zpWrkS3ngj7YokKf8ZtiRV6Npr4cADYdKktCuRpPzmMKKkCp13HrRqBXvumXYlkpTf7NmSVKHGjWHQoOT55Mkwbly69UhSvrJnS9JGxZhMmJ8zJ7lbsaioat/f574+AIwaPCrTpUlSXjBsSdqoEOCBB2D58qoHLUmSYUtSJeyxx5rnf/879O8PP/jBxr+nrEdr9Bej12rbwyWp0DhnS1KlffttsjTE9denXYkk5Q97tiRV2o47wuuvw157bfrcUYNHwbBh9PnmP7B0KaPunQrXXJPtEiUp59izJalKOnWCevVg/vykl2vu3A2cOGwYnHUWLF2atL/4ImkPG1ZjtUpSLnAjakmbZdQoOOIIeOIJOPTQCk5o1SoJWOtq2RKmTs1ucZJUwza2EbXDiJI2S58+MGUK7LRT0o5x7c2t4xfTCBV947RpNVCdJOUOhxElbbayoHXGGckdiwsXJu0Y4YLGdzOUK9f/phYtaq5AScoBhi1J1RIjfP110sv1m9+sCloXwC3zz2BO3WasNVGhYUMnyUsqOIYtSdUSAvzrX3DOOXD33VCnDtxyC5x/Ptx0bxNCy5bJSS1bwp13rtkDSJIKhBPkJWVEjEnQKrNkCdSvn149klSTNjZB3p4tSdVWNnRY3m9/C6Wl658nSYXGsCWpWlbP0Vo1dFhamjzeeit07gwLFqx93tChqZYrSTXOsCWpWkKAJk1WzdG6KWnfeCN06QLjxyeBq3wgmzPHHi5JhcV1tiRV29Cha6+zVacOjBsHAwbA88+vmctVPpBJUqGwZ0tSRqwboOrUgeeeW/vYlVcatCQVnoyErRDCYSGESSGEySGESyp4/cAQwrgQwooQwnGZuKak3FbRpPk+fRxClFR4qh22QghFwO3A4UB74KQQQvt1TpsGDAYerO71JOW+iibNH3ccfPBBctzAJamQZKJnqycwOcY4Jca4DHgYOKr8CTHGqTHGD4DSit5AUu1S0aT5Rx9N2o0bJwugTpmSdpWSVDMyMUF+F+DLcu3pwD4ZeF9JeWzdSfMhJMHr009h//2he3fYbbdUS5RUABYtSnYKS1NOTZAPIZwVQigJIZTMmDEj7XIkVdO6k+FDgLZtYdKkZPNqgOXLa74uSYVhwQLYay949dV068hE2PoK2LVcu/mqY1UWY7wzxlgcYyxu1qxZBkqTlIu22y55nDgxN/4ilFQ7NWoEf/kLtGqVbh2ZCFtjgTYhhNYhhC2AE4ERGXhfSbXc1ltD69awyy5pVyKpNrn/fnj55eT5wIHQokW69VQ7bMUYVwDnAC8AE4FHY4zjQwhXhRAGAoQQeoQQpgPHA/8bQhhf3etKyn+77AIvvgi7756033473Xok5b9ly+CGG+D229OuZI0Qc/Qe7OLi4lhSUpJ2GZJqyFNPwdFHw9NPw5FHpl2NpHxTWprclFNUBN99B9tsA/Xr19z1QwjvxBiLK3otpybISypcAwYk/xM9/PC0K5GUb0pL4dRT4dxzk8C1ww41G7Q2xbAlKSfUrQu//nXyv9I5c+DCC5NbtiVpU+rUSaYlNG+ediUVcyNqSTln1Cj461+TVef33z/taiTlqq++gqVLkzX7rrsu7Wo2zLAlKeccfTR89hn84AdJe/Fi2HLLVEuSlGNKS5P5nUVFMHZsbm9y7zCipJxUFrReeSX5X+u776Zbj6TcUqcO/M//wN1353bQAnu2JOW4Fi1g332T9bgk6a67ksczz4T99ku3lsqyZ0tSTtt9d3jiiWRj65Ur4cEHk7uNJBWe0tLk74Onn86vvwfs2ZKUNx59FAYNgu23h0MOSbsaSTVl4cLkP1tbbw2PPZYs65DrQ4flGbYk5Y0TT0yCVv/+STvG/PoLV1LVlZYmf+abNIFnn4Wttkq7oqpzGFFS3ghhTdD6/HPYZx/46KN0a5KUXXXqJPOzfv3r/P3PlT1bkvLS3LnJoqe5tEq0pMz5+9+hVSvo3RtOPz3taqrHni1JealrV/jgA2jTJmn/61/JcIOk/Ld0KVx7bbK4cW1g2JKUt+qs+hvsjTfgiCPW3BIuKT/NmQMrViQ91i+9BA88kHZFmWHYkpT39tsvuUPpjDOSdj7dEi4pMWcO7L03XH550v7BD6BevVRLyhjDlqS8F0Kyj2K9eskt4gccACNGpF2VpKpo0gROOSXZrqu2MWxJqlXmz096tho0SLsSSZuybBn89rcwZUrS/v3vkx0jahvDlqRaZaed4PXX1yx6+txz8H//l25Nkir29dfJXMtnnkm7kuwybEmqdcrW4pk/PxmW+O//TrceSWsbPz55bNUKPv4Yzjsv1XKyzrAlqdZq3BhGjYKbb07aCxe6PISUthdegE6d4KmnkvaOO6ZbT01wUVNJtVqnTsljjEkv17JlySa2+boStZSvyrbX6tsX/vjHNbtBFAJ7tiQVjMMPh0MPNWhJNe2ll+Cgg2DBguSu4YsvhoYN066q5tizJakghJDsr1bm9ddh5Ei44gqo69+EUtZ9/z3MnAmNGqVdSc2zZ0tSQXr6afjHP2Dx4rQrkWqniRPhkUeS5/36wbhxyYT4QmTYklSQrr0WSkqSSfQrVybhy5XnpcwZMiS5E3jJkqRdVJRuPWkybEkqWNtumzwOGwYDByZ3LkrafN9+mwwVAvzP/8Dbb7vAMDhnS5IYNAi23hr69EnaM2ZAs2apliTlncWLoXt36NULHnrIP0Pl2bMlqeAVFSX7sYWQrDbfrh3ccEPaVUn5Ydmy5HHLLeG665KbTrQ2w5YkldOkSXLX4pFHJm3ncUkb9tFH0KYNjBmTtAcNgvbt060pFzmMKEnlNGiQLLhY5oILksB1882uzyWtq1Ur6NABttoq7Upymz1bkrQBMUKdOknIMmhJiWeegWOPTe7ibdQo2ex9773Triq3GbYkaQNCgBtvhJtuStoTJ8LgwTBrVqplSamaMwc+/zy5kUSVY9iSpE0o69V6661k1fmVK9OtR6pJpaVw110wfHjSHjQIxo6FnXZKt658YtiSpEoaPBgmT4YddkiGGK+8Ej7+OO2qpOyKEe6+Gx59NGmH4BZXVWXYkqQqKNs8d9o0uOWWZL6KVNssXAh/+EPyWFSU/D4v23pHVWc2laTN0LIlfPIJNG2atEePhnnz4Ec/SrcuKRPefz/puW3XDo47DrbbLu2K8pthS5I20w47rHl+000wYQIcdhjUq5deTdLmev99+OADOOUU2H9/+PRT2H33tKuqHRxGlKQMeOyxZKilXj1Yvhyuvhq+/z7tqqTKu/ZauOwyWLo0aRu0MsewJUkZUK8e7LFH8nzMmGQIpmxVbSkXLV2a9MhOm5a0b7456d2qXz/VsmolhxElKcP69k3uUiwLXw89lOwbd9RRLo6q3PHtt3DppclSJv/937DjjmlXVHsZtiQpC9q0SR5jhL/+Nen5OuqodGuSxoyBUaOSzaJbtIDx4x0urAkOI0pSFoUAL7+c9G6FkMzjGjQoWa+rsvrc14c+9/XJWo0qHM8/nyxQOm9e0jZo1QzDliRlWd26a+5cfPddePZZWLQoaceYXl2q/b7+Go45Bv7zn6R96aUwaRJsvXW6dRUahxElqQb17QtffQVbbZW0L7oIZs6Ee+9dfz5XWW/W6C9Gr9UeNXhUzRSrvLV0aTLRfeutk6HCKVNg332TjaNV8wxbklTDyoJW2fOlS9cErW++gZ13Tqcu1Q6/+U2yd+FrryXhauLEZBV4pcdhRElK0e9/D7fdljyfMiWZtHzvvUm7z9RRdH1vFL1b9qZ3y968clrSHjo0tXKVo8aOTTaMBujWLelBXbEiaRu00mfYkqQcse22yaKShx6azOX67LNk/8VPP01ev+CCpD1njnO9tMZLL0HPnvDEE0n7tNOSfQ3dySB3OIwoSTli221Zq9eqQQPYYgv4+v+N4mtgNHD++clClK7XVbhWroT//V/Yfns44QTo0wfuuCPZKkq5yZ4tScpRf/kLvP762sdatzZoFar585PHOnXgvvvgySeTdlER/OIXa88FVG4xbElSjtpiC/jHP9Y+9vjjyRBiaSm88sqaeTqq3a69NlkTa8mSJGy/8AIMG5Z2Vaosw5Yk5aAY18zROv/8JFSdf36yAvgFF8DIkckk6KeeSrtSZcPMmUnA+vbbpH3AAUnv1bJlSXvbbe3hzCfO2ZKkHBQCNGmy9hytm25KXmvSJJmn8+CDMGBAcuyee+CZZ+D++11LKV8tWwYLFkDTpjB7drIA6U47weDBSdg64IC0K9TmykjYCiEcBtwCFAF3xxivXef1+sD9QHdgFvCTGOPUTFxbkmqroUOTHq6yHoyywFXWPumkNecuXgxz566Zt/P009C8ebIMgHLfihXJfppHHJHspdm2LUybBrvumnZlyoRqDyOGEIqA24HDgfbASSGE9uuc9jPg+xjjHsBNwHXVva6kaho2DFq1SmbbtmrlBJAcFR5c+3MKD1b8OZ19drIEQAhJQDv//GQNrzLTprlcRK655RY444zked268NvfwnHHrXndoFV7ZGLOVk9gcoxxSoxxGfAwsO7e9kcBf1/1fDjQLwRHm6XUDBsGZ50FX3yR/Av8xRdJ28CVWzbzcwoBSkrghhuS9pw5sMceyRwgSN7K4FXzJk5M1r8q+9nPmQMzZiRLOUASmPv1S608ZVEmwtYuwJfl2tNXHavwnBjjCmAusF0Gri1pc1x22ZqdkMssWpQcV+6oxufUtGkSsCDpNbn1VvjRj5L2++8nK9W/8UaG69ValiyB559PhncB3n4brroKPvkkaQ8Zkgz3usJ77ZdTdyOGEM4KIZSEEEpmzJiRdjlS7TVtWtWOKx0Z+pwaNYJf/hI6dlxzrLg4WUoAYMSIJIjNnLmZdQpIeqwmTEg2Ggd47z04/HB49tmkfdxx8N13sOeeSdvxncKRibD1FVB+ZLn5qmMVnhNCqAtsQzJRfi0xxjtjjMUxxuJmzZploDRJFWrRomrHlY4sfU5duyZbu+y4Y9L+/vtkA+xtt03at90GZ57pUOOmxAiTJsHHHyft779PAu3ddyftHj3gX/+CY45J2lttteZnrMKSibA1FmgTQmgdQtgCOBEYsc45I4DTVj0/Dng5Rv8YS6m55hpo2HDtYw0bJseVO2roczrttGSOV9lw1syZyfSwsp6Xc86B885bc34hL6T6n//Ayy+vaffpA1dfnTxv2hSGD4ef/SxpFxUlW+hsuWWNl6kcU+2lH2KMK0II5wAvkCz9cE+McXwI4SqgJMY4Avgb8EAIYTIwmySQSUrLoEHJ42WXJUNSLVok/4CXHVduSOlzKn8XIyQbGtcp91/z4mLYf/9kOyGADz5IbpjceuuslpWKf/4zCZ4XXJC0L7kkWWbjrbeSMPrAA8mvvcyxx6ZSpnJcyNUOpuLi4lhSUpJ2GZKkcmKEK69M1oQ65ZSkl6tsTtiNNyavX3opHHUU7Ldf0l65Mpmkn4u++SYZBjzooKR9/fVJ79RbbyXts85KtkX65JMkXE2cmAwF7rRTejUrN4UQ3okxFlf0Wk5NkJck5bYQkjvqTjklaZeWwsMPw6mnJu2ZM5M7H995Z027QYNkhXtIVka/6KJk8jgkN1e+996aTZarK0aYN2/NcgpTp8JDD63Z5uaf/4RevZI7BQHuuivZ9mjx4qS9ww7JBPYVK5L2rbeuCVoA7doZtFR1hi1J0marWxcGDkwm3QM0a5YEpzPPTNohJENvXbok7a+/ToYfp05N2hMnJqvcjxqVtN9+O9mOqHy7a9c14e3115OJ5x99lLRfeCEZxiubpP7gg7DNNjBlStIeNQpOPhmmT19Tc1HRmnA3aBC8+uqanrfTTku2PCprN2jgXYOqvhzt2JUk5auiojWT7bfffs0Eckju1lu0aM2djrvtltwZ2aNH0t5226SX7Ac/SNpbbAEtW665T6B+/aT3qSwMNWsGvXsnoQiS4Hb99UlggyQITpiQbF0EyZyq8vOqdt99zRIYUrY4Z0uSJKmanLMlSZKUEsOWJElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZZNiSJEnKIsOWJElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZZNiSJEnKIsOWJElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZZNiSJEnKIsOWJElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZZNiSJEnKIsOWJElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZZNiSJEnKIsOWJElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZVK2wFUJoGkL4dwjh01WP227gvOdDCHNCCM9U53qSJEn5pro9W5cAL8UY2wAvrWpX5HrglGpeS5IkKe9UN2wdBfx91fO/A0dXdFKM8SVgfjWvJUmSlHeqG7Z2jDF+s+r5/wE7VvP9JEmSapW6mzohhPAisFMFL11WvhFjjCGEWJ1iQghnAWcBtGjRojpvJUmSlBM2GbZijAdv6LUQwrchhJ1jjN+EEHYGvqtOMTHGO4E7AYqLi6sV3CRJknJBdYcRRwCnrXp+GvBUNd9PkiSpVqlu2LoW6B9C+BQ4eFWbEEJxCOHuspNCCGOAx4B+IYTpIYRDq3ldSZKkvLDJYcSNiTHOAvpVcLwE+Hm5dq/qXEeSJClfuYK8JElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEkFrM99fehzX5+0y5CkWs2wJUmSlEXVWvpBUn4q680a/cXotdqjBo9KpyBJqsXs2ZIkScoie7akAlTWg2WPliRlnz1bkiRJWWTPllTA7NGSpOyzZ0uSJCmLDFuSJElZZNiSJEnKIsOWJElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZZNiSJEnKIsOWJElSFhm2JEmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZZNiSJEnKIsOWJElSFhm2JEmSsijEGNOuoUIhhBnAF2nXkWe2B2amXYTW4meSe/xMcpOfS+7xM6maljHGZhW9kLNhS1UXQiiJMRanXYfW8DPJPX4mucnPJff4mWSOw4iSJElZZNiSJEnKIsNW7XJn2gVoPX4mucfPJDf5ueQeP5MMcc6WJElSFtmzJUmSlEWGrTwTQjgshDAphDA5hHDJBs45IYQwIYQwPoTwYE3XWIg29bmEEFqEEF4JIbwbQvgghHBEGnUWkhDCPSGE70IIH23g9RBCuHXVZ/ZBCGHvmq6x0FTiMxm06rP4MITwRgihS03XWIg29bmUO69HCGFFCOG4mqqttjBs5ZEQQhFwO3A40B44KYTQfp1z2gCXAj+MMXYAflPTdRaaynwuwOXAozHGbsCJwP/UbJUF6T7gsI28fjjQZtXXWcBfa6CmQncfG/9MPgd6xxg7AX/AOUM15T42/rmU/T13HTCyJgqqbQxb+aUnMDnGOCXGuAx4GDhqnXPOBG6PMX4PEGP8roZrLESV+VwisPWq59sAX9dgfQUpxvgqMHsjpxwF3B8T/wGahBB2rpnqCtOmPpMY4xtlf3cB/wGa10hhBa4Sf1YAzgUeB/w3ZTMYtvLLLsCX5drTVx0rry3QNoTwegjhPyGEjf5vRRlRmc9lKPDTEMJ04DmSv7iUrsp8bkrPz4B/pV2EIISwC3AM9v5uNsNW7VOXZFikD3AScFcIoUmaBQlIPov7YozNgSOAB0II/vmTKhBCOIgkbF2cdi0C4Gbg4hhjadqF5Ku6aRegKvkK2LVcu/mqY+VNB96KMS4HPg8hfEISvsbWTIkFqTKfy89YNScixvhmCKEByb5jdsmnpzKfm2pYCKEzcDdweIxxVtr1CIBi4OEQAiR/bx0RQlgRY3wy1aryiP+zzi9jgTYhhNYhhC1IJlqPWOecJ0l6tQghbE8yrDilBmssRJX5XKYB/QBCCO2ABsCMGq1S6xoBnLrqrsR9gbkxxm/SLqqQhRBaAP8ETokxfpJ2PUrEGFvHGFvFGFsBw4FfG7Sqxp6tPBJjXBFCOAd4ASgC7okxjg8hXAWUxBhHrHrtkBDCBGAlcJH/O8yuSn4u/0UypHsByWT5wdEVhbMqhPAQyX88tl81V+5KoB5AjPEOkrlzRwCTgUXA6elUWjgq8ZkMAbYD/mdVL8oKN0LOvkp8LqomV5CXJEnKIocRJUmSssiwJUmSlEWGLUmSpCwybEmSJGWRYUuSJCmLDFuSJElZZNiSJEnKIsOWJElSFv1//xehYAK3o2YAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "poly = pd.DataFrame(index=X.index)\n",
    "\n",
    "poly[0] = 1\n",
    "poly[1] = np.cos(X)\n",
    "poly[2] = np.sin(X)\n",
    "\n",
    "# technically no need to recreate poly, however we do since the indices may change, or may increase or decrease\n",
    "\n",
    "model = sm.OLS(Y, poly)\n",
    "res = model.fit()\n",
    "coef = res.params\n",
    "\n",
    "continuation = (poly * coef).sum(axis=1)\n",
    "\n",
    "exercise = (k - df[N-1][ITM])\n",
    "\n",
    "x = np.linspace(.5, 1.5, 100)\n",
    "\n",
    "y = 1 * coef[0] + np.cos(x) * coef[1] + np.sin(x) * coef[2]\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "\n",
    "plt.plot(x,y, linestyle=\":\", color=\"blue\")\n",
    "\n",
    "plt.scatter(X, Y, color=\"red\", label=\"Y (discounted exercise later)\")\n",
    "\n",
    "plt.scatter(X, continuation, label=\"continuation\", marker=\"x\", color=\"blue\")\n",
    "plt.scatter(X, exercise, label=\"exercise now\", marker=\"+\", color=\"green\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we are done: we know for each option when the best time to cash it in was. So we know the **stopping rule**. With that, we can compute the option value, by discounting the cash flow by exercising the option at the optimal point."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 64-bit",
   "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}