FiMa I - Homework Assignment No. 5 Solution

FIMA1_WS1920_sheet5_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following exercises will guide you towards implementations of pricing tools for \n",
    "European options in the CRR-model.\n",
    "\n",
    "Concerning the return of those exercises, please send an email with your jupyter (ipynb) file to either\n",
    "\n",
    "          [email protected]   or   [email protected]\n",
    "          \n",
    "until Wednesday, 20 November 2019, 12:15 p.m. Indicate your group members in the email."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Setting** We consider the Cox-Ross-Rubinstein model with time horizon $T>0$.\n",
    "Using the notation from the lecture, we assume\n",
    "$$\n",
    "\t-1 < d < r < u.\n",
    "$$\n",
    "to make the model arbitrage-free and complete. We want to tackle the following tasks:\n",
    "\n",
    "**Setting for Exercise 5.2.1**\n",
    "    \n",
    "Let $T=3$. Calculate the *Price*\n",
    "$$\n",
    "    P^{call}(0):=S^0(0)\\mathbb{E}^{\\mathbb{Q}}\\left[\\frac{(S^1(T)-K)_+}{(1+r)^T}\\right]\n",
    "$$\n",
    "at time $t=0$ for the plain-vanilla call option\n",
    "for parameters \n",
    "$$\n",
    "    S^0(0)=1,\\,\\, S^1(0)=120,\\,\\, d=-0.5,\\,\\, u=0.5,\\,\\,  r=0.1,\\,\\,  K=120\\,\n",
    "$$\n",
    "with unique martingale measure $\\mathbb{Q}$.\n",
    "\n",
    "**Setting for  Exercise 5.2.2**\n",
    "\n",
    "Now we consider a so called *up-and-out-Call* with *knock-out* price $B$ ($B=200$).\n",
    "This claim is a European Call on the asset $S^1$ with maturity $T$ and strike $K$ if\n",
    "the asset price never exceeds the knock-out price $B$ within the time period $[0,T]$.\n",
    "If, however, the asset price is larger than $B$ at any time point $t\\in \\{0,\\ldots,T\\}$,\n",
    "the holder gets nothing.\n",
    "\n",
    "a) Give a mathematical description of the discounted payoff of the\n",
    "up-and-out call-option (in the notebook (see cell below)). \n",
    "\n",
    "b) Calculate the discounted price $P^{k.o.}(0)$.\n",
    " Use the same parameters as in 5.2.1 with $T=3$. Discuss your findings also for $T=2$ and $T=30$. \n",
    "\n",
    "     \n",
    "**The goal of the programming exercise is to write Python code which does this computations for general $T$.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Toy example\n",
    "\n",
    "For your convinience, let us give the following python code that calculates (not particularly efficiently) $C(0)$ from Exercise 5.2.1 for $T=1$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#In this cell we use the \"Code\"-environment to code something in python. \n",
    "#More precisley we just load all necessary packages.\n",
    "import numpy as np \n",
    "#This package contains a lot of methods to implement mathoperations\n",
    "#in a fast way\n",
    "import matplotlib.pyplot as plt \n",
    "#This package is used for generating plots\n",
    "from scipy.stats import bernoulli \n",
    "#This package is fast for simulating bernoulli random variables\n",
    "from scipy.special import comb\n",
    "#This package allows to compute n over k fast\n",
    "import time \n",
    "#Used to stop the compution time\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#input variables:\n",
    "S0 = 120\n",
    "a = -0.5\n",
    "b = 0.5\n",
    "r = 0.1\n",
    "B = 200\n",
    "K = 120\n",
    "T = 1\n",
    "#Unique equivalent martingale measure:\n",
    "p_star = (r-a)/(b-a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the price via the backwards formula from 5.1 iv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculated value is 32.73.\n",
      "The time needed was 0.00000 seconds.\n"
     ]
    }
   ],
   "source": [
    "#Backward formula of the expected value for T=1 only!!!\n",
    "tic_1 = time.time()\n",
    "pi_1 = p_star*max(S0*(1+b)-K,0)/(1+r)+(1-p_star)*max(S0*(1+a)-K,0)/(1+r)\n",
    "toc_1 = time.time()\n",
    "print( \"The calculated value is %.2f.\" %(pi_1))\n",
    "print(\"The time needed was %.5f seconds.\" %(toc_1-tic_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the price via Monte Carlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The following function provides independent samples of the call at time T\n",
    "def sample_call_value_one(a,b,S0,K,r,M):\n",
    "    random_values = bernoulli.rvs(size=(M,1),p=p_star)\n",
    "    up_moves = random_values*(1+b)\n",
    "    down_moves = -1*(random_values-1)*(1+a)\n",
    "    moves = up_moves+down_moves\n",
    "    return ((np.maximum(S0*np.prod(moves, axis=1)-K,0))/(1+r)**1)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The last calculated value is 32.71.\n",
      "The time needed for the last run was 0.04000 seconds.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8ZGWV8PHfqb2yp5NeQnqnG5Cl2ZpNEATZnNdBEBQQEBDBBbdhBoVhRgXH133fAFGZmVfEBdEWxVYBQUGBbqAXlqZXeqXTW/ak1vP+cW8l1emqSipVqaok5/v51Ofe+9ztPJVKnbrPfe69oqoYY4wxo+UpdwDGGGPGN0skxhhjCmKJxBhjTEEskRhjjCmIJRJjjDEFsURijDGmIJZIjCkSEXlRRN5c7jiMKTWx60iMGRsicq+qXlPuOIwZa3ZEYkwRichsEfm6iFS500eJyJfLHZcxY8lX7gCMmShEZBPwPuAB4H+BxUAv8Gl3/r1ADzAXOB14CXi3qq4vfbTGFI8dkRgzthJAMm36cuB2oBFYB3yuHEEZU0yWSIwproOAi4GrgMeAu4Fb0ub/SlWfUdU48BPgmNKHaExxWdOWMcW1XVX/B0BEUNVVwM1p819PG+8FakoZnDFjwY5IjBkj1mPLTBaWSIwxxhTEEokxxpiC2AWJxhhjCmJHJMYYYwpiicQYY0xBLJEYY4wpiCUSY4wxBZkwFyQ2Nzfr3Llzyx2GMcaMK8uXL9+tqlML2caESSRz585l2bJl5Q7DGGPGFRF5rdBtWNOWMcaYglgiATQaRROJcodhjDHj0qRPJNHNm1l33vl0/enP5Q7FGGPGpUmfSPytrXiCQXbffRd2lb8xxuRv0icS8Xppuv56Ii+9TM8TT5Q7HGOMGXcmfSIBqL/gn/Ed1MLuO+2oxBhj8mWJBBC/n6brrqPv+efpffbZcodjjDHjiiUSV8PFF+NtbmbPXXeXOxRjjBlXLJG4PKEQTddcTc+TT9K3alW5wzHGmHHDEkmahssuw1NXx+677ip3KMYYM25YIknjralhypVX0v3nR4isXVvucIwxZlwYUSIRkTkicrY7HhaR2rENq3war7oSqapi990/KHcoxhgzLgybSETkeuCXQKq9Zybw67EMqpx8jY00Xnopnb/7HdHNm8sdjjHGVLyRHJHcCJwKdAKo6lpg2lgGVW5Trr0G8fnYc88Pyx2KMcZUvJEkkoiqRlMTIuIDJvRVe/5p06i/+B10PPggsZ07yx2OMcZUtJEkksdF5N+BsIicA/wC+O3YhlV+TdddhyaT7P7Od+1qd2OMyWEkieQWYBewCng/8HvgP8YyqEoQmDmTxiveTfsvfsHOz/1fNJksd0jGGFORhn1CoqomgR+4r0ll+i23IF4fe3/8Y+J7dnPQF7+IJxAod1jGGFNRhk0kIrKRDOdEVHX+mERUQcTjYfonP4GvuZm2L3+ZLfvamfmdb+OtqSl3aMYYUzFG8sz2xWnjIeCdwJSxCady/Pr5bXx56Rq2t/dxUMNcPvuhT9Jy91d57T3vYfZdd+GbOrXcIRpjTEUY9hyJqu5Je21T1W8AZ5UgtrL59fPbuPVXq9jW3ocC29r7uHHXDLb86x1EN25i07uvIPraa+UO0xhjKsJILkg8Lu21WEQ+AAx7ZbuIhETkGRFZISIvisjtbvkP3bKVIvJLEcnYTiQit4rIOhFZIyLn5V2zAnx56Rr6Yvs/w70vluD2nfXMuffHJLu62HTZ5ez7xS/sWe/GmElvJL22vpr2+jxwPPCuEawXAc5S1aOBY4DzReRk4F9U9WhVXQRsBj48dEURORy4DDgCOB/4noh4R7DPotje3pe1PHz00cz56X0E5szh9f/8FBsvvIjuv/6tVKEZY0zFGUnT1plpr3NU9XpVXTOC9VRVu91Jv/tSVe0EEBEBwmS+uPHtwP2qGlHVjcA64MQR1qlgBzWEc5YH581jzk/vo/Ub3yDZ38+W669n8/uup3/Nq6UK0RhjKkbWk+0iclOuFVX1a8Nt3D2KWA4sAL6rqk+75T8G/gl4CfjXDKu2Av9Im97qlg3d/g3ADQCzZ88eLpwRu/m8Q7n1V6v2a94K+73cfN6h6fum7vzzqDnrTPbddx+7v38nGy+6iIaL30HTddcRmDu3aPEYY0wly3VEUjvMa1iqmlDVY3Bu9HiiiBzpll8LHAS8DFyaYVXJtLkM279bVRer6uKpRexFdeGxrXz+HUfR2hBGgNaGMJ9/x1FceOwBuQxPIEDTNdewYOkfmHLVlbT/+jesP/+tzjmU+39GoqOjaHEZY0wlklLd/kNEPg30qOpX0srOAG5W1bcNWfZWAFX9vDu9FPiMqv492/YXL16sy5YtG5PYR+rXz2/j7gef5vAXn+L8bc/R2r4DCQSoOess6t9+ATWnnYb4/WWN0Rhj0onIclVdPPyS2Y3kgsQQcB3Oie9QqlxV3zvMelOBmKq2i0gYOBv4kogsUNV17jmSfwZeybD6EuA+EfkazpHLQuCZEdapLFJdhvtiQV5aeCa/XPBmjujZwW2+jXie/gtdf/gDnro6qk85herTTqXmtNPwt7SUO2xjjCnYSC5I/F+cL/vzgDuAK3CapIbTAvy3e57EA/wc+B3wVxGpw2m+WgF8EEBELgAWq+qnVPVFEfk5zjmUOHCjqlZ0P9sDugyL8GLNQdzccDB/e+J2uv/6N7oefYSevz1J19KlAAQOPpia006l+tRTCS9ahLehoUzRG2PM6A3btCUiz6vqsSKyUlUXiYgfWKqqFXVRYrmbtubd8ruM3c8E2PiF/zMwrapE16+n+69/Y+PDfyawegX+ZByA6IxWmhcfS/ioIwkdtYjQ4W/AEwpl2KoxxhRHSZq2gJg7bHdPlr8OzC1kpxPRQQ1htmW4/mRoV2IRIbhgAQ93hbl1wwySsy/msH2vcci+LbyhYyuLn3oa/0MPOQt7vQTmzSW4YCHBhQsILlxIcMFCArNnIb6R/OmMMWbsjeTb6G4RacS5dfwSoAb4zzGNahwaSZfhdANNYb4AK6YuZMXUhYDTQ+wv7z2S/tWr6Vu1isira+l/6SWnOcw9epRAgMC8eQTmziUwZ44znDuXwLy5eBsacE4/GWNMaeS6jmS6qu5U1XvcoieACX/H39FKdQ0evNFjmJvPOzRjl2HIffW8f9o0/GedRe1Zg62Hyb4+Ius3EFm7lsjatUQ3bCCyZg1djzwC8fjAcp66OvwHHeS8WlqcYasz7psxA19Tkx3NGGOKKtc3ygoRWQX8FHhAVe2CiGFceGxr1sQx1EibwlI84TDhI48gfOQR+5VrPE5s2zaimzY5r9c2E9u+ndi2bfQ++yzJrq4hG/Lga2rCN306vmnT8E2fhm/qVHzNzc6rqQlvUzO+5iY7P2OMGZFciaQVp8vuZcDnReTvOElliapm/jltRizfprBsxOdzmrfmzIEzzjhgfqKri9j2HcS2byO+cyfxtjZiO3cS39lGbMsW+pYvz3rRpKe6Gu+UKXinNOJrnIK3sdEZnzIFb0MD3vp659XQgCc1tAd/GTPpjOiCRBEJAG/FSSpnAo+o6hVjHFteyt1razT2f+ZJ7qawsZSMREjs2UN8zx7iu3c747ud6cS+fST27iW+b9/AuEajWbcl4TDe2lq89XV4auvw1tXhqavFW1uHp7YGb00NnprawfHaWjzVNXhrqt3xasRbsvtzGjPplarXFqoaFZGXcK4fOR44vJCdGkc+TWFjyRMM4nHPqwxHVUn29JLsaCfR0UGi3R2mxts7SHR3kezodI6G2naSXLeORGcnye5uSCaH3YdUVTlJprraeVVVDb6qq5Bw2BkPp8qdaQmH3bIwnlDInR4cF89IbnZtjMlXzkQiIrNx7oV1OVAN3A+8XVVHckGimYBEBG9NNd6aavyt+SVBVUX7+kh0dZPs6SbZ1eWMd7vT3d0kurp5dcPrrF67HXp7aIjFWRDvZ0p3N8neXufV00Oyr29ESWm/2IPBwQQTDDrDUAgJhZzpUAhPKIgEQ0goiCcYctcJIoGgW+bODwbc8SASCCCBIJ5gwBlPlQWDiN9vvejMhJer19ZTOOdJfgHcoKrjq93IVBwRQdwjC5iWcZlfP7+NW3euou+I/c8dDb1ppqqikQjJvj6SPb1on5tk+vpJ9vWi/f37j/f2oRG3rL8P7Y+Q7O9H+/pIRiLEu7pIRvqdcneo/f1oLJYpzPzq7fe7yWbIa6DcjycQQPxDElFqvt/vjvudZVLTfj/i87nlQ8rccXypct/+89whqXE7WjMFyHVEcivwhJbqro4TWKWcCxkPsj2d8stL1+z3nomIewQRgsbGMYtHk0k0EnGSVmrY349GY2g0bToSdaajUWe5aMxZL+qWx2Iko1F3OuYMY+4wGiXZ20cy2n7g/NjgOGP5NE6Px0koPt9gckl74fchPrfc680+7fUiPi/4fIjX54x7U9vZv1x8Pmee1ws+74HlPi94vc42vd6BsvRxPN60srR5Xo8znUqS6UOv16mvHSkWTdZEoqqPlzKQiWrwZo7Ol8C29j5u/dUqgKImk4mSrHJdX1MO4vEg4TCEw5SzC8Cvn9/GbQ+sIBaJ4ksm8CUT1HiVfz93AecsnDKQdIjHnfHUMBZzEtNA+ZDlYnFnXtwpW7u9nec37CbSH6HWLyyaUc2s2oC7TGq5hDMejZLs64XUNhIJZ7vuuCbig8smEhCLOcM8myTHSkI8JMWD1+fFl0paHicB4fU4iW2/JORFPPsvNzBMreNJH2ZYbmB5j5ME04brdvfx9GvtdEYTVIUCnLJwKocdVO9sWzzOsuIZ3L7HM7C+r7mJ2re8pWzvpV2ZNsZG+gu7EKVKVqWQ7/U1E1GmHwVfXrqGnriC10/U6zyKoBP4/HMdvO3c44uyD8D5HB2Wu1mxUJpMDiac1HjCSTiDiSgJiVR5wh1POoksmUwrSwwmsIF1kk4SSySdbSYTA/Ne3LKPP6/eQTIex6tJPJokIHD6gikc3FQFiQSaTDjrJhNOIky6240PmZdIDmxbE3GIxUgOmT8Qx9DyZAKS6u4vSSwWpyoa5y1uTB5VPCuV3SN8T8NHH22JZCIrxS/sUiSrUinW9TXjVbYfBUP/vinDfY5yJowh+wj5PaP6HOV7NCweD7jniUrt1i88yraFB75nDzWEefKW8t2H9tQvPHrgDyhVZtaH+OvNZ7jJLIEmFZLOUV16mXjLe45rTB+1a0rzC7vSmoMKke+tZlLy/TKr1KbAbD8KvCIkMpyuzPU5ypaUsiWM0SSrSj4azvQ3rtT/lYz7F2FbZ8TpHFH6kPKS64gk9TjdQ4ETcG7YCM7DqJ4Yy6AmklL8wp5ozUH5Xl+T75fZWHz5FSsxZftCS6gS9nvz+hxlS0rZEkY2qc9Rtia3Yh4NF+t9zPY3bqjys6/3wJ545f5fGe//w7lOtt8OICJ/BI5T1S53+jM4XYLNCIz2F3Y+StkcVM5f8tn2nevLDA5870f75Zdt/7kSU6b959pHti+U1rTYM+0/U3m+v7Ibwn4i8WTGz1E5m9wgd4LPJ8EFfZ68E/Jo5fq8DC0f7026I3mw1SvA0aoacaeDwApVPawE8Y3YeLxFSjHl+oIfq195MDYnY/Pd97/87IWMDxVLLTN0nWxffqmHkI3kSy59/19euibjl3+2L+Zc71e+73Gu5UcTF2ROfBnb8CFrk1trjnMO2WIO+T0ZjxaGS6KZtpXrb/z1S48Z8x9D2eK6+PhWHli+La/3PrW9fH5E5KMYt0gZSSK5DXgX8CCgwEXAz1X1/xay42IbbSL59fPbiD70CWbH1hPwepg9pYrmmuAYRFgeu7sjbNjdQzLt7+wRYX5zNQCb9/YSTSRHVPfnNu8jmjiw62Zq3Xy2la9c+wYyzstXKu5M75dHIJ488H8l4PXkve+A18Nxs7Nf+7K7OzLi93K4v0mx/vb/2Lgn6zyPSF77yBZzLtn2kdr+SA333ucr298q3zrmiivb//DUmgC7uqMkVXkpOYc74u8Z1Q+7ktxrS1U/JyIPA29yi65V1ecL2WmlSP1quFnj4HG+jDbs7gGYMMlk897e/T6AAElVNu3pIakMzBtJ3bP9Y6TWzbStVAyFJphc+14wtSbjP9rQeqfL9MWUSoaZ3q8MOWRg//kmk+GWba4Jjvg9yvW+pLaR7f3P5++QrY7ZfkQAWT8To0n6mf4mwyWRbH/j0ciUMKB4dcy1fLbP5M6uyAHLlqu35ki7/1YBnar6YxGZKiLzVHXjWAZWCql21Dt4z37lrT1hnvxIeR9JX6zmqIuzPEs+m1x1/9c8mzca9mZoQmnz8vnT8m8Ky7bv1gYn3r9laSfP53zDice2ckK+75e7br5NNU9eW5zPV8735dqzaAaai7Cfzdma0C44iuOObT1gH6d+4VG2RTLE1ROGKvJqcsvVTJXvOaUTi9Wk2+b+fSMZ/r456pizKTDLZyLf/+Fy9EAbNpGIyKeBxTi9t34M+IH/B5w6tqGNvUrtCljMXkXZ/tGy2d7elzWJZfvCzPaP3t534D/ZaH8xDXcyMltPr2zrZFs+2/uV60R0tg4VufZfLKU6SVvMJ4B+/dJjMsb8mQuOyLiPbD8IDsqRxHP9jfOVb++3XHXMdo4k198r22dyNF3Cx8pIjkguAo4FngNQ1e0iUpt7lfGhUrvcFbNLZb6/luvD/mGT2Ej/0bMZTaIeTe+30ayT7f3K9iWX2lauL62xPLFbil6B6fsqxhNAh4s53x8EubZVDPl+Xoer4+I5U4rymRxNUhorIznZ/oyqnigiz6nqcSJSDfxdVReVJsSRGc3J9nL2QsplXpZD2VSvonzl0wspZ3NMEXvhlPMq4uFU6sWK40mx/7fK+TfJ1mNtNL3yRmsi9Nr6N2AhcA7weeC9wE9V9VuF7LjYCum1VWlXRGf74Bb7CzhTXbJ1pR0uieXbZda+mCe+8ZSQVZWEJkhqcr9hIpngdyu38bnfv0h/PA6SAJSQX/j4OQejmuTepzayq7uX5toA7z5pJqcumHLAdoYO99tfMvNyqVem6fT1k5qkpbqFKw+/clR1L0kicXd0DnAuzvfJUlX9UyE7HQuluI6kVEcw5TxSKnYSG09fJqVUqe+LqhJPxoklY8Q1TjwZJ5FMEE864zGN7Ted0ISzrDueWj61XKosNX9gXXW3O3SYvh33izJ93YQmBrfjrpNxesi6yWTSmZ82nvoST2gCzet0dnkJgle8zkPmxItHPBzVfBT3nHfP6LZXoiOSL6rqJ4crK7dSJJJSHSlA+b5oKrW5b7xLapJoIko0GWXJC5v53MOr6I9HEUmAJ07Qn+T602dz4vy6geViidjAl3o0ESWWjA2+ErGBL/BU2cB4IrbfsgOJIW2YbTyhY/jMkyy84sXn8eEVL16PF5/48Hq8B5SnvjR9Ht9+ywwM09b3iGe/cp/HLRsy7vV4B8Y94mHN6z389dU9dPYlqAsFOfsNMzh2VhMej2dgWwPLezz7TQ/MH1KePkxf18OB28m1TipxFPs5KqVKJM+p6nFDylZOhHMk+Sr2uYtKVam/loshlozRH+8fePUl+ojEI/Qn+okkIvuN98f7iSai9Cf2H0YSESKJyH7jqelUWSwRI5J0ymIJ59d9Mfk9fvwePz6Pzxn3+vGJj4A3MFjmzs82nW3ewEt8B5SlvogHtuF+oQ9dZ+BLPb0sLTmkvsx94quYB0xN1h9RY3pBooh8EPgQcLCIrEybVQs8VchOx6tK7eVVbMXqNjmcbAkrqUl6Y730xHrojfc6r1gvffE+emOD073xwbK+eB+98V4nOcT7BpJEX6yP/sRg4hjtF7pXvAS9QYLeIAFvgJAvRMAbIOhxpqt91TQGG50yd5mAJ+AM3XG/10/AE+CzD70KSR+qXlAfqBdVH6I+fvmB0wl4A/g9/v3WS33Zp77wK+XLdyKZSI9jKLVc3X/vAx7GOcF+S1p5l6ruHdOoKtR4v7FaKSQ1SU+sh85oJ13RLjoj7jDaSXes23lFu3nx9TaefW07ydp+wvX97PP28x/LI/zX6hiRZO+I9+cTH2F/mLAvTJWvirDPGa8J1DDVN5WwL0zIFyLkDe03HvI5r6A3SNgXHkgSqbKQN0TQFxwo93mK9+ieu377KNs6MjeRHjPtmKLtx+SnUq8rGw9y3f23A+gQkW8Ce9Pu/lsrIiep6tOlCrJSlLLPfjalbHaKJ+O0R9rZ27+X9v529kb2sq9/H+397bRH2umIdtAeaacz0klHxBnvinYNe+Iy7AvTH/GjgSCSDKKJEBqvJZkModEa3n/KYdQGaqnyV1Htq6bKX0WVr2q/YSpx+N2nBY4n9oOkMk2WFoexMJJzJM/j3EZe3WkPsGzoeZNymwx3/y1GG66q0hXrYlfvLnb27mRX7y529e1iT98e59U/OGyPtGfdTm2glvpAPfXBehqCDdQF66gP1FMXrKMuMPiqDdTu96ryV+H3+CfN+aZsJvJ5qPHKzpGM3kiO10XTso2qJkXEHtFbBiNpw40kIuzo3sGOHue1vXv7wHhbbxttvW30xQ/81VXtr6Yp1ERTuIl59fNYPGMxTaEmGkONNIYamRKaQkOwgcZQIw3BhoKbeib7r79SnYcyI1cJLQ7j1Ui+DTaIyEeB77vTHwI2jF1IJpvBttoEEtiLJ7AbT2A3ewK7ed/S+9jUuYmdvTv3W8cjHqaGp9JS3cJhUw7j9JmnM71qOlPDU5lWNY3pVdNprmom7CvtF7g175hKZAl+dEaSSD4AfAv4D5znkTwC3DCWQRlHUpNs69rG2va1rGtfR8PcJ4nIdjzBXc71BwMLVtEXX8CJM05kVt0sWmtaaaluoaW6henV0/F7Ku88gv36M2biGNGV7aPasEgI59nuQZyE9UtV/bSI/ATnbsIx4Bng/ap6wA2ZxPmmTD2vdLOqXpBrf+P9HImqsr1nO6t2r2L1rtWs3rOal/a8tF8zVIN/OnvbpxDrm0YyOo1ktJlgcjqfv/Bk+wI2xoxKSc6RiMghOM1a01X1SBFZBFygqv81zKoR4CxV7RYRP/A39wFZPwFSN4W5D3gfg81m6fpUdcL2hVRV1rWv46ntT/H0jqdZvXs1+yL7AAh4AhzWdBgXLbiIw6YcxsENB3Nww8FU+6sHT9J22K94Y0xlGEnT1g+Am4G7AFR1pYjcB+RMJO4J+m530u++VFV/n1pGRJ4BZo4i7nGpI9LB33f8nae2PcWT25+krbcNgHn18zhj1hkc1XwURzQfwSENh2Tt1mptuMaYSjOSRFKlqs8MuZJ2RJcHi4gXWA4sAL6bfu2Je5RyFfCxLKuHRGSZu68vqOqvM2z/BtzzNbNnzx5JSCXXHe3mT6/9iSXrl/Bc23MkNUltoJZTWk7h1NZTeeNBb2RG9Yxyh2mMMaM2kkSyW0QOxjnRjohcAuwYycZVNQEcIyINwIMicqSqrnZnfw94QlX/mmX12e5DtOYDj4rIKlVdP2T7dwN3g3OOZCQxlUIimeAfO/7BkvVLeHTzo/Qn+plTN4frj7qe01pP48jmI4t6pbQxxpTTSL7NbsT5sj5MRLYBG4Er8tmJqraLyF+A84HV7uN7pwLvz7HOdne4wV33WGB9tuUrQUekg3tfvJcl65bQ1tdGbaCWCw6+gAsWXMCi5kV2fyRjzIQ0bCJR1Q3A2e6TET2pW6UMR0SmAjE3iYSBs4Evisj7gPOAt6hqMsu6jUCvqkZEpBnn+fBfGlmVSi+aiPLTV37KXSvvoifWw5ta38QnD/4kZ8w6g6A3WO7wjDFmTI2k11YT8GngNEBF5G/AHaq6Z5hVW4D/ds+TeICfq+pDIhIHXgP+7v5C/5Wq3iEii4EPqOr7gDcAd4lI0l33C6r60ijrOGZUlT++9ke+sfwbbO3eyqkHncpNi2/ikMZDyh2aMcaUzEiatu7HuR7kYnf6CuBnOEcYWanqSpzmqKHlGfepqstwugKjqk8BR40gtrJ5oe0FvrLsK6zYtYKFjQu58+w7ObX11HKHZYwxJTeSRDJFVT+bNv1fInLhWAVU6VSVbz//bX6w6gc0h5v5zCmf4cIFF+L1eMsdmjHGlMVIEsljInIZ8HN3+hLgd2MXUuWKJ+Pc8fc7eHDdg1y04CJuOfEWqvxV5Q7LGGPKaiSJ5P3ATcD/utNeoEdEbsK5wLBurIKrJH3xPj7x+Cf4y9a/8P5F7+fGY260XljGGMPIem3VliKQStYR6eAjj36EF9pe4LaTbuOywy4rd0jGGFMxPMMtICLXDZn2uteBTAqv97zO1Q9fzerdq/nKGV+xJGKMMUMMm0iAt4jI70WkRUSOAv4BTIqjlA3tG7jq4avY2buTO8++k3PnnlvukIwxpuKMpGnr3SJyKc4t3XuBy1X1yTGPrMySmuSmv9xELBHj3vPv5dAp9sAlY4zJZCRNWwtxbqz4ALAJuEpEJnxXpUc3P8r6jvV88sRPWhIxxpgcRtK09VvgP1X1/cAZwFrg2TGNqsxUlbtX3s3s2tmcO8eas4wxJpeRJJITVfURcPr6qupXgQl9QeKT25/k5b0vc91R19mFhsYYM4ysiUREPgGgqp0i8s4hs68d06jK7Acrf8CM6hn88/x/LncoxhhT8XIdkaT3c711yLzzxyCWirB853Kea3uOa464JutTCo0xxgzKlUgky3im6QnjByt/wJTQFN6x8B3lDsUYY8aFXIlEs4xnmp4QXtzzIk9uf5KrDr+KsC9c7nCMMWZcyHUdydEi0olz9BF2x3GnQ2MeWRncs/Ieav21XHropeUOxRhjxo2siURVJ1V3pfXt6/nz5j9zw6IbqA1Migv3jTGmKEbS/XdS+OGqHxL2hbnyDVeWOxRjjBlXLJEAW7q28PuNv+eSQy6hMdRY7nCMMWZcsUQC3Lv6Xjzi4Zojril3KMYYM+5M+kTS1tvGg+se5MIFFzKtalq5wzHGmHFnJE9InNBq/DV89NiP8pY5byl3KMYYMy5N+kRS5a/imiOvKXcYxhgzbk36pi1jjDGFsURijDGmIKI6Me52IiK7gNeGWawZ2F2CcCrVZK7/ZK47TO76W91zm6OqUwvZyYRJJCOVFYDFAAAf3UlEQVQhIstUdXG54yiXyVz/yVx3mNz1t7qPfd2tacsYY0xBLJEYY4wpyGRLJHeXO4Aym8z1n8x1h8ldf6v7GJtU50iMMcYU32Q7IjHGGFNklkiMMcYUZNIkEhE5X0TWiMg6Ebml3PHkQ0R+JCJtIrI6rWyKiPxJRNa6w0a3XETkW249V4rIcWnrXO0uv1ZErk4rP15EVrnrfEtEJNc+SklEZonIYyLysoi8KCIfm2T1D4nIMyKywq3/7W75PBF52o3tZyIScMuD7vQ6d/7ctG3d6pavEZHz0soz/m9k20epiYhXRJ4XkYdyxTXR6i4im9zP5Qsisswtq8zPvapO+BfgBdYD84EAsAI4vNxx5RH/6cBxwOq0si8Bt7jjtwBfdMf/CXgY55HIJwNPu+VTgA3usNEdb3TnPQOc4q7zMPDWXPsocd1bgOPc8VrgVeDwSVR/AWrccT/wtFuvnwOXueV3Ah90xz8E3OmOXwb8zB0/3P3cB4F57v+DN9f/RrZ9lOE9uAm4D3goV1wTre7AJqB5SFlFfu5L/qEo0wfxFGBp2vStwK3ljivPOsxl/0SyBmhxx1uANe74XcDlQ5cDLgfuSiu/yy1rAV5JKx9YLts+yvw+/AY4ZzLWH6gCngNOwrla2eeWD3y+gaXAKe64z11Ohn7mU8tl+99w18m4jxLXeSbwCHAW8FCuuCZg3TdxYCKpyM/9qJu2RMQjIu8a7fol1gpsSZve6paNZ9NVdQeAO0w9TCVbXXOVb81QnmsfZeE2VRyL86t80tTfbdp5AWgD/oTzK7pdVePuIukxD9TTnd8BNJH/+9KUYx+l9A3gE0DSnc4V10SruwJ/FJHlInKDW1aRn/tRJxJVTQIfHu36JSYZyiZqv+dsdc23vKKISA3wAPBxVe3MtWiGsnFdf1VNqOoxOL/OTwTekGkxd1is+pf9fRGRtwFtqro8vTjDohOu7q5TVfU44K3AjSJyeo5ly1rHQk+2/0lE/k2cE6JTUq8CtzkWtgKz0qZnAtvLFEux7BSRFgB32OaWZ6trrvKZGcpz7aOkRMSPk0R+oqq/Gia2CVf/FFVtB/6C0wbeICKp5wmlxzxQT3d+PbCX/N+X3Tn2USqnAheIyCbgfpzmrW/kiGsi1R1V3e4O24AHcX5EVOTnvtBE8l7gRuAJYLn7WlbgNsfCs8BCtydGAOdE3JIyx1SoJUCqB8bVOOcOUuXvcXtxnAx0uIenS4FzRaTR7YVxLk677w6gS0ROdnttvGfItjLto2TcmH4IvKyqX0ubNVnqP1VEGtzxMHA28DLwGHBJhtjSY74EeFSdxu4lwGVuz6Z5wEKck60Z/zfcdbLtoyRU9VZVnamqc924HlXVK3LENWHqLiLVIlKbGsf5vK6mUj/3pT6BVK4XTq+GV3Hal28rdzx5xv5TYAcQw/klcR1OO+4jwFp3OMVdVoDvuvVcBSxO2857gXXu69q08sXuh3Q98B0G73iQcR8lrvtpOIfcK4EX3Nc/TaL6LwKed+u/GviUWz4f58twHfALIOiWh9zpde78+Wnbus2t4xrcHjq5/jey7aNM/wNvZrDX1oSvu7v/Fe7rxVRslfq5L+gWKW6TwwdxuqeCc9h9l6rGRr1RY4wx40qhieQenL7t/+0WXQUkVPV9RYjNGGPMOFBoIlmhqkcPV2aMMWbiKvRke0JEDk5NiMh8IFHgNo0xxowjvuEXyelm4DER2YBzsmcOcG3BUY1Cc3Ozzp07txy7NsaYcWv58uW7tcBnto86kYiIB+jD6Up3KE4ieUVVI4UENFpz585l2bJK7HlsjDGVS0ReK3Qbo04kqpoUka+q6ik4XRONMcZMQoWeI/mjiFycuv3weNQfS/D4q7vYsre33KEYY8y4VGgiuQnnYp2IiHSKSJeI5LoPUsXp7I9x9Y+eYemLr5c7FGOMGZcKOUciwBGqurmI8ZTctNoQLfUhVm7tKHcoxhgzLhVy91/FuZHYuHdUaz0rt7aXOwxjjBmXCm3a+oeInFCUSMro6FkNbNrTS0ev3dnFGGPyVWgiORP4u4isd58TvEpExl0PrqNa6wFYtc2at4wxJl+FXpD41qJEUWaLZjqJZOW2dk5b2FzmaIwxZnwZ1RGJiJwFoKqvAR5VfS31Ao4vZoCl0FAVYE5TFSu32BGJMcbka7RNW19JG39gyLz/GOU2y+qo1npr2jLGmFEYbSKRLOOZpseFo2c2sK29j93dZbnDizHGjFujTSSaZTzT9LhwVOo8iXUDNsaYvIz2ZPt8EVmCc/SRGsednleUyErsyNZ6RGDl1g7OOmx6ucMxxphxY7SJ5O1p418ZMm/o9LhQE/Rx8NQau8LdGGPyNKpEoqqPFzuQSrBoZj1PvLrbeZj9+L0PpTHGlFShFyROKIta69ndHWFHR3+5QzHGmHHDEkmaRbMaAKx5yxhj8lCURCIi1cXYTrkd3lKHzyPWc8sYY/JQUCIRkTeKyEvAy+700SLyvaJEVgYhv5dDptfahYnGGJOHQo9Ivg6cB+wBUNUVwOmFBlVOR8+qZ+XWDpy75BtjjBlOwU1bqrplSFGi0G2W01GtDXT0xdhsj941xpgRKTSRbBGRNwIqIgER+TfcZq7xKnUn4BV2wt0YY0ak0ETyAeBGoBXYChzjTo9bh86oJeDzsHKLnXA3xpiRKOSZ7V7gKlW9oojxlJ3f6+HwljpW2gl3Y4wZkUKe2Z5g/1ulTBiLZtazelsHiaSdcDfGmOEU2rT1pIh8R0TeJCLHpV5FiayMFs1soDeaYMOu7nKHYowxFa/QR+2+0R3ekVamwFm5VhKR84FvAl7gHlX9wpD5NwHvA+LALuC97tMXS+LotBPuC6fXlmq3xhgzLhWUSFT1zHzXcc+tfBc4B+cE/bMiskRVX0pb7Hlgsar2isgHgS8BlxYSaz7mT62hKuBl1dZ2Ljl+Zql2a4wx41KhRySIyP8BjgBCqTJVvSP7GpwIrFPVDe769+OcaxlIJKr6WNry/wCuLDTOfHg9wpGt9dYF2BhjRqDQW6TciXOk8BGch1q9E5gzzGqtQPpFjFvdsmyuAx7Osv8bRGSZiCzbtWvXiOMeiUWt9by0o5NYIlnU7RpjzERT6Mn2N6rqe4B9qno7cAowa5h1Mj3oI2P3KBG5ElgMfDnTfFW9W1UXq+riqVOn5hH28BbNaiAaT7Lm9a6ibtcYYyaaQhNJnzvsFZGDgBjDP2p3K/snm5nA9qELicjZwG3ABaoaKTDOvC1qdU642w0cjTEmt0ITyUMi0oBzxPAcsAm4f5h1ngUWisg8EQkAlwFL0hcQkWOBu3CSSFuBMY7KnKYq6kI+u6W8McYMo9BeW591Rx8QkYeAkKrm/AmvqnER+TCwFKf7749U9UURuQNYpqpLcBJTDfAL95G3m1X1gkJizZeIsHjuFB59pY1IPEHQ5y3l7o0xZtwoKJGIyHsylKGq/5NrPVX9PfD7IWWfShs/u5C4iuXqN87l6h89w6+f38alJ8wudzjGGFORCm3aOiHt9SbgM0BJjxzG0ukLmznioDrufHyD3S7FGGOyKCiRqOpH0l7XA8cCgeKEVn4iwo1nLmDj7h7+sPr1codjjDEVqSjPbE/TCyws8jbL6rwjZjC/uZrvPrbOnppojDEZFHpB4m9FZIn7eghYA/ymOKFVBq9H+MAZB/PSjk4ef7W4Fz0aY8xEUOgtUr6SNh4HXlPVrQVus+JceGwrX//zq3zvL+t586HTyh2OMcZUlELPkTye9npyIiYRgIDPw/Vvms8zG/eybNPecodjjDEVpdCmrS4R6czw6hKRzmIFWQkuO3EWjVV+vveX9eUOxRhjKkqhJ9u/DtyCc9PFmcAngf9S1VpVrSs0uEpSFfBx7anzePSVNl7eMaFypDHGFKTQRHKeqn5PVbtUtVNVvw9cXIzAKtHVp8ylOuDl+3ZUYowxAwpNJAkRuUJEvCLiEZErgEQxAqtE9VV+rjx5Dg+t3M5re3rKHY4xxlSEQhPJu4F3ATuBNpznkby70KAq2XWnzcPn9XDn4xvKHYoxxlSEQnttbVLVt6tqs/u6UFU3FSm2ijStLsQlx8/kgeVbWbvTnlVijDGjSiQicr2ILHTHRUR+JCIdIrJSRI4rboiV5+NvWUhd2M8H/t9yuiPxcodjjDFlNdojko/hPHsE4HLgaGA+cBPwzcLDqmzT6kJ8+/Jj2bi7h1seWGm3TjHGTGqjTSRxVY25428D/kdV96jqn4Hq4oRW2U45uIl/O+9QHlq5g/9+alO5wzHGmLIZbSJJikiLiISAtwB/TpsXLjys8eEDpx/M2W+Yxud+/zLPbd5X7nCMMaYsRptIPgUsw2neWqKqLwKIyBnApOnO5PEIX33nMcyoD3HjT55jT3fJHy1vjDFlN6pEoqoPAXOAN7jPIUlZBlxajMDGi/oqP9+/4nj29ET5+M9esAdgGWMmnVF3/1XVuKruG1LWo6rdhYc1vhzZWs8dFxzBX9fu5puPrC13OMYYU1LFfrDVpHXpCbN45/Ez+faja1myYnu5wzHGmJKxRFIkIsJnLzySxXMa+ehPn+fOx9dbt2BjzKRQcCIRkVYReaOInJ56FSOw8Sjk9/K/153E2xa18IWHX+HfH1xNLJEsd1jGGDOmCnpCooh8Eefk+ksM3qxRgScKjGvcCvm9fOuyY5nTVMV3H1vP1n29fO+K46gN+csdmjHGjIlCH7V7IXCoqlq/1zQej3DzeYcxe0oVtz24mnfe+Xd+dM0JHNQwaS6xMcZMIoU2bW0A7Kd2FpeeMJt7rz2Rbfv6uPC7T7Jya3u5QzLGmKIrNJH0Ai+IyF0i8q3UqxiBTRSnLWzmgQ+9Eb/Xwzu+9xRf+9OrROIT9pEtxphJqNCmrSXuy+RwyPRaHvrIaXz2oZf41iNreXjVDr54ySKOm91Y7tCMMaZgMlG6qC5evFiXLVtW7jCG9diaNm771Sp2dPbz3lPn8a/nHkJVoNB8bowxoyMiy1V1cSHbKKhpS0QWisgvReQlEdmQehWyzYnuzEOnsfRfTufKk+bww79t5LxvPMHf1u4ud1jGGDNqhZ4j+THwfSAOnAn8D/C/hQY10dWG/Hz2wiP52Q0n4/N4uPKHT3PlPU/z9IY95Q7NGGPyVmgiCavqIzhNZK+p6meAs4ZbSUTOF5E1IrJORG7JMP90EXlOROIickmBMVask+Y38fDH3sStbz2MV17v5NK7/8Gld/2dJ9fttqvijTHjRqGJpF9EPMBaEfmwiFwETMu1goh4ge8CbwUOBy4XkcOHLLYZuAa4r8D4Kl7I7+X9ZxzMXz9xFp962+Fs2tPDFfc8zcXff4rH1rRZQjHGVLxCE8nHgSrgo8DxwJXA1cOscyKwTlU3qGoUuB94e/oCqrpJVVcCk+b+IuGAl/eeNo/Hbz6Tz154JDs7I1z742c5+2uPc89fN7CvJ1ruEI0xJqOCugup6rMAIqKqeu0IV2sFtqRNbwVOGs3+ReQG4AaA2bNnj2YTFSfk93LVyXO4dPEsfvPCNu57ZjP/9buX+dIf1vDWo2Zw+YmzOWneFESk3KEaYwxQ+L22TgF+CNQAs0XkaOD9qvqhXKtlKBtV+42q3g3cDU7339Fso1IFfB7euXgW71w8i5d3dHL/M5v51fPb+M0L25k/tZp3LZ7F2xa1MLOxqtyhGmMmuUKbtr4BnAfsAVDVFcBwd//dCsxKm54J2AM8cnhDSx23v/1Invn3s/nKO4+msSrAFx5+hdO++BgXf/8p7n1yI21d/eUO0xgzSRV8JZyqbhnSzDLc/T+eBRaKyDxgG3AZ8O5C45gMwgEvlxw/k0uOn8nmPb38duV2frtiO5/57Uvc/tBLnDyvibcd3cKZh06zG0QaY0qm0ESyRUTeCKiIBHBOur+cawVVjYvIh4GlgBf4kaq+KCJ3AMtUdYmInAA8CDQC/ywit6vqEQXGOqHMbqrixjMXcOOZC1i7s4vfrtzBQyu2c9uDqwFYOK2GMw6ZyhmHTuWEuVMI+b1ljtgYM1EVdIsUEWkGvgmcjXPu44/Ax1S15FfWjZdbpIwlVeXVnd088eounli7i6c37CWaSBLyezh5fhNvPLiJE+c1ceRBdfi89nBMY0xxbpFi99qawPqiCf6xcQ+Pr3ESy4ZdPQBUBbwcP6eRk+ZN4cR5TSyaWW9HLMZMUsVIJKNq2hruVvGq+tHRhWOKKRzwcuah0zjzUOca0baufp7duI9nNu7h6Y17+cofXwXA7xXe0FLH0TMbOHpWA8fMqmd+cw0ej3UxNsYMb1RHJCISBVYDP8fpcbXfN46q/ndRosuDHZHkr703yjMb9/L8lnZe2NzOqm0ddEfiANQGfRzZWs8RB9VxRGsdh7fUc/DUamsSM2aCKdsRCdACvBPnee1x4GfAA6q6r5BgTGk1VAU494gZnHvEDAASSWXDrm5e2NLOiq3trNzawf/84zWicecGAwGfh8Nm1HJ4Sx2HzqjlkOm1LJxew9SaoF0gacwkVvA5EhFpBS4HbgI+qaplufuvHZGMjXgiyfpdPby0o4MXt3Xy0o5OXtzeSUdfbGCZ+rCfQ6bXsHB6LQum1jBvajUHN9fQ2hjGa81jxlS0ch6RpAI4DieJnAM8DCwvZHum8vi8Hg6dUcuhM2q56FinTFXZ1RVhbVs3r+7sYm1bN2t3dvG7lTv2SzABr4fZTVXMa65mfnM1s5uqmNVYxewpVRzUECbgs2YyYyaC0Z5svx14G841I/cDt6pqvJiBmcolIkyrCzGtLsSpC5oHylWVvT1RNuzuYeOuHme4u5sNu3p4fM0uoonBe3B6BFrqw8yeUsWsKamh+2qsorkmYM1lxowToz3ZngQ2AH1uUWojAqiqLipOeCNnTVuVLZlUdnb1s3lPL5v39rJlby9b9vXx2p4etuzrY1dXZL/lw34vrY1hWupDtNSHmFEfdofudF2I+rDfko0xBSpn09a8QnZqJh+PR2ipD9NSH+ak+U0HzO+LJti6r5ct+3rZvMdJMlv39fJ6Z4RXd+6irSvC0N88QZ+H6XVOUpleH2JGXZBptSGm1QWZWhNkaq0zXRf2WcIxZgyNKpGo6mvFDsRMbuGAl4XTa1k4vTbj/FgiSVtXhNc7+tjR0c/rHf3s7OxnZ2eE1zv7Wbm1nT929BOJH/gIm4DXw9TaIM01AZprgs6rdnC8qTrAlJoAU6oDTKkKWBdnY/JU8E0bjSkFv9dDa0OY1hw3o1RVOvvj7OqKsKsrQltX/8D4rq4Iu7ojbO/oZ+W2Dvb2REkkMzfr1of9TnKpDtBYHaCxyk+jm2ScaaesoSpAQ5WfhrDfko+Z1CyRmAlDRKgP+6kP+1kwrSbnssmksq83yu7uKHt6IuztibKvJ8qenih73eG+nihb9vaycmuUfT2x/ToLDFUb9NFQ7ey7IRxw4qjyD8RTH/ZTF/JTG/JRF/ZTNzD0W+81M+4V2v33Y6r6zeHKjKk0Ho/QVBOkqSYIZG5OS6eq9EQT7OuJsq83yr7eGO29Udp7Y7T3xtjXG6W9N0pHX4yOvhjbO/rodMdjidwdWkJ+z0CiqQ/7qXMTT23IR03QR62bgFKvmqCfmqAzrybkozroJeize6WZ8in0iORqnLv/prsmQ5kx45qIDHx5z5oy8qdSqiq90QSd/TE6++J09sfoShvv6I0NzOvoc8Z3dvaztq2Lrv44Xf3xrE1w6QJez0BSqQ44Cac6lWyCznh1wOsMg4PLVQd9VLnlVQGnrCroJeD1WAcFM2KjvY7kcpyHUc0TkSVps2pxn5ZojHESUOrLu6U+//VVlb5Ygu7+OJ39cbr6Y3RH4vREnCTTE4nTHYnTFYnTPTCdoCcSZ093lM17eumKxOmNxOmJDvfMuUE+jxAOeKkKeKkK+Aj7vVQHvYQDPqr8Tnlqfjjgc5fzEvJ7Cade7nRV4MBpv51TmlBGe0TyFLADaAa+mlbeBawsNChjjENEqAr4qAr4mFZX2LaSSScppZJPTyRBbzRObzRBTzROb8QdRp1leqOD8/uiCXqjCTp6o+xwx/tizvz+WPZzR9n4PELY7yUU8BLyewaST9AdhvweQn4vId/geDBV7nMS0sAyblnQ7yXoGywLpq3r84gdYY2hQrr/vgacIiLTgRPcWS/bFe7GVCaPZ/DoaFoRt5tKUH0xJ+H0p40PlMUT9EWTbuJJJaEE/bEkkdTy7rLtvVH6Y0n64862+mNJ+mOJjF27R1x3YSCxBH1egn4PQZ+HgM9DwOsOfU6TXtDnvvZb1jtY7nOmA760bbjbCfq9A9sLps33u2UTNaEVerL9ncBXgL/gXNX+bRG5WVV/WYTYjDHjQHqCGkvJpBJNJPdLLk6ySZU5yaY/liASSxJx56UPB+a7w1hCicaTRONJOvpi7rgzPxJ3klxqvBhEGExc3v2T0AHjXicB+X0e/F4ZmB5ITF4ZmN9SH+Ltx7QWJcbRKPQv/x/ACaraBiAiU4E/A5ZIjDFF5fEIIY+3LE/zVFUi8STRhJN0IvHU0Elag+WJ/eanyqNDplPbig1dJuHM647EiSWSxOJKLDG4TCyRdJKfO51y7OyGcZ1IPKkk4toD2Fk0Y8yEIiLuuZfK6WatqsSTTqIZQce+MVVoIvmDiCwFfupOXwr8vsBtGmOMGYaI4Hebt8qtoESiqjeLyDuA03DOkdytqg8WJTJjjDHjQjHOjj0JxHBuJf9MEbZnjDFmHCnoUbsi8i7gywz22noTUJZeWyKyC6dLci7NwO4ShFOpJnP9J3PdYXLX3+qe2xxVnVrITgpNJCuAc4b22lLVowsJaqyIyLJCH+Aynk3m+k/musPkrr/VfezrXuhZGuu1ZYwxk9xY9Np6uMBtGmOMGUcmW6+tu8sdQJlN5vpP5rrD5K6/1X2MFXSO5ICNiXiBy1T1J0XbqDHGmIo2qvMZIlInIreKyHdE5FxxfBjYALyruCEaY4ypZKM6IhGR3wD7gL8DbwEagQDwMVV9oagRGmOMqWyqmvcLWJU27sVJKrWj2VapXsD5wBpgHXBLuePJM/YfAW3A6rSyKcCfgLXusNEtF+Bbbj1XAselrXO1u/xa4Oq08uOBVe4632LwB0bGfZS47rOAx4CXgRdxfqxMpvqHcC70XeHW/3a3fB7wtBvbz4CAWx50p9e58+embetWt3wNcN5w/xvZ9lGG98ALPA88NJnqDmxyP5cvAMsq+XM/2go+l2u60l7uB3E9MB/nyGkFcHi548oj/tOB49g/kXwp9cEHbgG+6I7/E07POQFOBp5O+3BscIeN7njqQ/gMcIq7zsPAW3Pto8R1b0n9U+A8gfNV4PBJVH8Batxxv/vldjLwc5zzkQB3Ah90xz8E3OmOXwb8zB0/3P3cB3G+JNe7/xdZ/zey7aMM78FNwH0MJpJJUXecRNI8pKwiP/ejrWAC6HRfXUA8bbyzHB+2YeI9BViaNn0rcGu548qzDnPZP5GsAVrc8RZgjTt+F3D50OWAy4G70srvcstagFfSygeWy7aPMr8PvwHOmYz1B6qA54CTcK5W9rnlA59vYClwijvuc5eToZ/51HLZ/jfcdTLuo8R1ngk8ApwFPJQrrglY900cmEgq8nM/qpPtqupV1Tr3VauqvrTxAh8IOiZagS1p01vdsvFsuqruAHCHqYfeZatrrvKtGcpz7aMsRGQucCzOr/JJU38R8YrICzjNm3/C+RXdroNPI02PeaCe7vwOoIn835emHPsopW8AnwBSD9/IFddEq7sCfxSR5SJyg1tWkZ/7sX2kWeXI9GzL4vV7rizZ6ppveUURkRrgAeDjqtqZ43GlE67+qpoAjhGRBuBB4A2ZFnOH+dYz04/JinhfRORtQJuqLheRN6eKMyw64eruOlVVt4vINOBPIvJKjmXL+rmfLLcz2Ypz0jZlJrC9TLEUy04RaQFwh6lb1WSra67ymRnKc+2jpETEj5NEfqKqvxomtglX/xRVbce5QerJQIOIpH4Ipsc8UE93fj2wl/zfl9059lEqpwIXiMgm4H6c5q1v5IhrItUdVd3uDttwfkCcSIV+7idLInkWWCgi80QkgHMibkmZYyrUEpzeGLjD36SVv8e9tudkoMM9PF0KnCsijSLSCJyL0+67A+gSkZPF+Zn/niHbyrSPknFj+iHwsqp+LW3WZKn/VPdIBBEJA2fj9GB7DLgkQ2zpMV8CPKpOY/cS4DIRCYrIPGAhzsnWjP8b7jrZ9lESqnqrqs5U1bluXI+q6hU54powdReRahGpTY3jfF5XU6mf+1KfQCrXC6dXw6s47cu3lTuePGP/KbAD57kvW4HrcNpxH8HpovcIMMVdVoDvuvVcBSxO2857cbr6rQOuTStf7H5I1wPfYbAbYMZ9lLjup+Eccq/E6Qb5gvu3nCz1X4TT9XWlG+On3PL5OF+G64BfAEG3POROr3Pnz0/b1m1uHdfg9tDJ9b+RbR9l+h94M4O9tiZ83d39r2Cw2/dtuT6T5f7cF/UWKcYYYyafydK0ZYwxZoxYIjHGGFMQSyTGGGMKYonEGGNMQSyRGGOMKYglEmOKSESecofiDj+TPm3MRGTdf40ZAyJyE86NTA8FosDjqvrH8kZlzNiwIxJjikhEugHUuQq/Gfgo8AdV/aOIvFlE/iIivxSRV0TkJ3akYiaCyXLTRmNKSkQ+jnPPpm8B54tICOfOBMcCR+Dc1+hJnPtJ/a1ccRpTDHZEYszY+Kaq3gP0qOptwJ/d8mdUdauqJnFu9zK3XAEaUyyWSIwZA+qefFTVz6RPA5G0xRJYq4CZACyRGGOMKYglEmOMMQWx7r/GGGMKYkckxhhjCmKJxBhjTEEskRhjjCmIJRJjjDEFsURijDGmIJZIjDHGFMQSiTHGmIL8fxIMcDDCW7HfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Value of pi via Monte Carlo\n",
    "n=10000\n",
    "mc_values = 50\n",
    "tic_2 = time.time()\n",
    "x_values = []\n",
    "pi_M = []\n",
    "pi_M_RMSE = []\n",
    "for i in range(1,mc_values,1):\n",
    "    M=n*i\n",
    "    if i == mc_values-1:\n",
    "        tic_2 = time.time()    \n",
    "    x_values.append([M])\n",
    "    values_call = sample_call_value_one(a,b,S0,K,r,M)\n",
    "    pi_M.append([np.mean(values_call)])\n",
    "    pi_M_RMSE.append([np.sqrt(np.var(values_call,ddof=1))/np.sqrt(n*i)])\n",
    "toc_2 = time.time()\n",
    "print( \"The last calculated value is %.2f.\" %(float(pi_M[len(pi_M)-1][0])))\n",
    "print(\"The time needed for the last run was %.5f seconds.\" %(toc_2-tic_2))  \n",
    "#The following will plot the calculated value for pi\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(x_values, pi_M, 'o')\n",
    "plt.plot(x_values, [pi_1]*len(x_values), '-')\n",
    "plt.plot(x_values, pi_1-1.96*np.array(pi_M_RMSE), '-')\n",
    "plt.plot(x_values, pi_1+1.96*np.array(pi_M_RMSE), '-')\n",
    "plt.title('i*n')\n",
    "plt.ylabel('Expected Value')\n",
    "#The following will plot the root mean square error\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(x_values, pi_M_RMSE, '-')\n",
    "plt.xlabel('i*n')\n",
    "plt.ylabel('Root Mean Square Error')\n",
    "plt.show()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 5.2.1 (4 points)\n",
    "\n",
    "Calculate the searched Price of Exercise 5.2.1 above for $T=3$ (rather than $T=1$ as assumed \n",
    "in the code above) via jupyther in three ways:\n",
    "\n",
    "    1. Use the backwards algorithm.\n",
    "    \n",
    "    2. Via a combinatorial argument.\n",
    "\n",
    "    3. Use Monte-Carlo Simulation, stop the time and plot your results as in the toy example.\n",
    "    \n",
    "    4. Why did we not ask you to implement explicit integration over $\\Omega=\\{-1,1\\}^T$ here? Compare in general the numerical complexity, when T=50, for all of the above methods.\n",
    "    \n",
    "    \n",
    "   \n",
    "Next you should calculate the call option price in the **$N$-step CRR model**\n",
    "($N \\in \\mathbb{N}$). Specifically, we want to investigate the\n",
    "convergence of the computed prices to the prices obtained by the\n",
    "corresponding limiting *Black-Scholes pricing formula* as the\n",
    "number of time steps $N$ becomes large. \n",
    "\n",
    "The $N$-step CRR model with time step $T/N$ for $T = 1$, riskless interest rate $r=0$ (for simplicity) and risky asset\n",
    "$S^{(N)}$ whose dynamics are described by\n",
    "\\begin{equation*}\n",
    "S^{(N)}_0 := (S^{(N)})^1(0):=s:= 100 > 0, \\quad S^{(N)}_k:=(S^{(N)})^1(k):= s \\prod_{i=1}^k (1 + R^{(N)}_i) \\quad (k=1,\\ldots,N),\n",
    "\\end{equation*}\n",
    "where $1+R^{(N)}_k$, $k=1,\\ldots,N$, takes the values\n",
    "\\begin{equation*}\n",
    "d_N := e^{-0.3 \\sqrt{T/N}} \\quad \\text{and} \\quad u_N :=\n",
    "e^{0.3 \\sqrt{T/N}}.\n",
    "\\end{equation*}\n",
    "Furthermore the martingale\n",
    "measure $\\mathbb{Q}$ with respect to the filtration\n",
    "$(\\sigma(R^{(N)}_1,\\ldots, R^{(N)}_n))_{n=1,\\ldots,N}$ is given by\n",
    "\\begin{equation*}\n",
    "  q_N :=\\mathbb{Q}[R^{(N)}_k = b_N] = \\frac{1-d_N}{u_N - d_N}, \\quad\n",
    "  1-q_N =\\mathbb{Q}[R^{(N)}_k = a_N] = \\frac{u_N - 1}{u_N - d_N} \\quad (k=1,\\ldots,N),\n",
    "\\end{equation*}\n",
    "where $R^{(N)}_k$ $(k=1,2,\\ldots)$ are i.i.d. under $\\mathbb{Q}$. Note that\n",
    "the (time discrete) price process of the risky asset\n",
    "$(S_k^{(N)})_{k=0,1\\ldots,N}$ can be discretized *in space* via\n",
    "the formula\n",
    "\\begin{equation}\n",
    "  S^{i,N}_k := s d_N^{k-i} u_N^{i}, \\quad 0 \\leq i \\leq k, \\quad 0 \\label{disc}\n",
    "  \\leq k \\leq N,\n",
    "\\end{equation}\n",
    "where $i$ denotes the number of upward movements until time step\n",
    "$k$.    \n",
    "\n",
    "5. Implement the backwards recursion algorithm from 5.1 iv) which allows you to compute the arbitrage-free price $P^{put,(N)}(0)$ at time t = 0 of a European put option with strike K = 90 in the above N-step CRR model with step sizes $$N = 50,60,70,...,1980,1990,2000.$$ Plot the deviation of the computed prices $(P^{put,(N)}(0))_{N =50,60,...,2000}$ from the price obtained by using the corresponding limiting Black-Scholes pricing formula for put options, i.e.\n",
    "$$\n",
    "    P^{Put}(s,K,r=0,T=1,\\sigma=0.3)=K\\mathrm{e}^{-rT}\\Phi(-\\mathrm{d}_-(s,K,r,T,\\sigma))-s\\Phi(-\\mathrm{d}_+(s,K,r,T,\\sigma))\n",
    "$$\n",
    "with\n",
    "$$\n",
    "    \\mathrm{d}_\\pm(s,K,r,T,\\sigma):=\\frac{\\mathrm{ln}\\left(\\frac{s}{K}\\right)+\\left(r\\pm\\frac12\\sigma^2\\right)T}{\\sigma\\sqrt{T-t}}.\n",
    "$$\n",
    "What do you observe?\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "For your convenience, see below the code lines that demand your attention."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backward formula for arbitrary $T$ (5.2.1 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculated value is 116.5178.\n",
      "The time needed was 0.000999689102 seconds.\n"
     ]
    }
   ],
   "source": [
    "#backward formula of the expected value for general T\n",
    "tic_1 = time.time()\n",
    "#start the time\n",
    "asset_val = S0*((1+b)**np.arange(T,-1,-1))*((1+a)**np.arange(0,T+1,1))\n",
    "#Vector which contains all possible outcomes of the asset at T\n",
    "call_val = np.maximum(asset_val-K,0)/(1+r)**T\n",
    "#Vector which contains all possible outcomes of the call at T\n",
    "#Using the recursion formula for pricing in the CRR model:\n",
    "for t in np.arange(T,0,-1):\n",
    "    call_val[0:t] = p_star*call_val[0:t]+(1-p_star)*call_val[1:t+1]\n",
    "pi_1 = call_val[0]\n",
    "toc_1 = time.time()\n",
    "#Ending the timer\n",
    "print( \"The calculated value is %.4f.\" %(pi_1))\n",
    "print(\"The time needed was %.12f seconds.\" %(toc_1-tic_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combinatorial argument for arbitrary $T$ (5.2.1 2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculated value is 116.5178.\n",
      "The time needed was 0.000000000000 seconds.\n"
     ]
    }
   ],
   "source": [
    "#Combinatorial argument from Exercise sheet 5.1 (vi)\n",
    "tic_comb = time.time()\n",
    "#start the time\n",
    "asset_val = S0*((1+b)**np.arange(T,-1,-1))*((1+a)**np.arange(0,T+1,1))\n",
    "pi_comb = 0\n",
    "for i in range(T,-1,-1):\n",
    "    pi_comb += comb(T, i)*p_star**i*(1-p_star)**(T-i)*(np.maximum(asset_val[T-i]-K,0)/(1+r)**T)\n",
    "toc_comb = time.time()\n",
    "#Ending the timer\n",
    "print( \"The calculated value is %.4f.\" %(pi_comb))\n",
    "print(\"The time needed was %.12f seconds.\" %(toc_comb-tic_comb))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo for arbitrary $T$ (5.2.1 3.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The following function provides independent samples of the call at time T\n",
    "def sample_call_value(T,a,b,S0,K,r,M):\n",
    "    random_values = bernoulli.rvs(size=(M,T),p=p_star)\n",
    "    up_moves = random_values*(1+b)\n",
    "    down_moves = -1*(random_values-1)*(1+a)\n",
    "    moves = up_moves+down_moves\n",
    "    return ((np.maximum(S0*np.prod(moves, axis=1)-K,0))/(1+r)**T)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The last calculated value is 115.67.\n",
      "The time needed for the last run was 0.72600 seconds.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4HNXZ8OHfsyutepcty3KRbclyxw1jXMG4AQZsMAabQBIIJZWEQAIp8KaQEAgB8iYQIMmXhNBxob/G2AZjY4wbLrhJluUiWZLVu7TlfH/sSkjy7mqbtCrnvq69tDuzM/PMaHaemXPOnBGlFJqmaZrmDUOwA9A0TdN6Hp08NE3TNK/p5KFpmqZ5TScPTdM0zWs6eWiapmle08lD0zRN85pOHprmIxH5UkQuCXYcmhYMou/z0LTAEJF/KaW+Eew4NK0r6CsPTfODiAwRkSdEJNLxebyIPBbsuDSts4UEOwBN66lEJA/4FrAaeAGYCtQBDznG/wuoBdKBOcAhYJVS6njXR6tpgaWvPDQtsKyArdXnlcCvgAQgB3g4GEFpWqDp5KFp/hkIXAfcDGwGngPubzV+jVLqc6WUBXgRmNj1IWpa4OliK03zT4FS6j8AIoJS6gBwX6vxha3e1wHRXRmcpnUWfeWhaQGiW1ppfYlOHpqmaZrXdPLQNE3TvKZvEtQ0TdO8pq88NE3TNK/p5KFpmqZ5TScPTdM0zWs6eWiapmle69E3CSYnJ6v09PRgh6Fpmtaj7N69u0Qp1c+fefTo5JGens6uXbuCHYamaVqPIiIn/Z2HLrbSNE3TvNZnk4dqagp2CJqmaT1Wn0we1R9+yLHZczAXFQc7FE3TtB6pTyaPsMxMbJWVVK5bF+xQNE3TeqQ+mTxMQ4cSeeGFVKxejbLZOp5A0zRNa6NPJg+A+OuXYz51irqdurWWpmmat/ps8ohZuBBDTAwVq98Idiiapmk9Tp9NHobwcGKXXEn1+g+wVlUFOxxN07Qepc8mD4D45ctRjY1UvvNOsEPRNE3rUfp08ogYO5aw0aOpfGN1sEPRNE3rUfp08gCIv+46Gg4douHQoWCHomma1mN4lDxEZKiIzHe8jxCRmM4Nq+vEXbUEMZmo0FcfmqZpHusweYjI7cAbwLOOQYOADu+uE5F/ikixiBx0Mu5eEVEikuz4LCLyZxHJEZH9IjLZu9XwnTEujpiFC6l85x1sDQ1dtVhN07QezZMrj+8CM4EqAKVUNtDfg+n+BSxuP1BEBgMLgFOtBl8OZDpedwDPeDD/gIlffh22qiqqN2zoysVqmqb1WJ4kj0alVEsvgiISAqiOJlJKbQHKnIx6AvhJu3lcA/xH2X0GxItIqgexBUTktGmEDh6si640TdM85Eny+FhEfgZEiMgC4HXgbV8WJiJXA/lKqX3tRqUBp1t9PuMY5mwed4jILhHZde7cOV/COH+eBgPx111L3Y4dNJ061fEEmqZpfZwnyeN+4BxwALgTeA/4hbcLEpFI4OfAg85GOxnm9OpGKfWcUmqqUmpqv35+PQirjbhly8BgoGLNmoDNU9M0rbfqMHkopWxKqeeVUtcrpZY73ndYbOXECGAYsE9E8rBXvO8RkQHYrzQGt/ruIKDAh2X4LDQlhejZs6lcsxZlsXTlojVN03ocT1pbnRCR3PYvbxeklDqglOqvlEpXSqVjTxiTlVKFwFvALY5WV9OBSqXUWW+X4a+45ddhKS6mZuvWrl60pmlaj+JJsdVU4ELHazbwZ+C/HU0kIi8D24EsETkjIre5+fp7QC6QAzwPfMeDuAIu5pJLMCYnU/Hqa8FYfK9S/vrrnLr9Dt3lvab1UiEdfUEpVdpu0JMishXndRetp1vZwfj0Vu8V9ibBQSWhocRfv5zSvz1L05l8TIOc1tlrHVA2G6XPPY/59GlqP/mE6Llzgx2SpmkB5kmx1eRWr6kichfQa+4wby9hxQoQoeLVV4MdSo9Vt3MX5tOnQYSy/74Y7HA0TesEnhRbPd7q9XtgCrCiM4MKptDUVGIum0fFG29ga2wMdjg9UuWa1Riio0n61reo/eQTGnNPBDskTdMCzJPWVpe2ei1QSt2ulDraFcEFS8LKlVjLy6lev77D79bv28eJ65ZTu317F0TW/Vmrq6la/wGxV15J4tdvQUJDKX/ppWCHpWlagLms8xCRe9xNqJT6U+DD6R4ip0/HlJ5O+YsvEXf11S6/t27XKSK//2MGl+eTe+u3KLntbi699w6/lq1sNsTQczs7rnr/fVRDA/HXXUtIcjKxV1xO5Zo19Pvh3Rijo4MdnqZpAeLuKBXTwavXEoOBhFUrqd+3j/ovv3T6nXV78/nkj39jcHk+f5q0gr39Mhnw9yfYdu+DKKvVp+XWbNtG9uw51O3Z40/4QVW5eg2mjBGEjx8PQMLXvoatro7KtR32palpWg8ivt3v1z1MnTpV7dq1q1Pmba2qInvOXOKuWkLqb35z3vjLH1rHr1f/imMJQ/j5jNsxKBt3HniTq098SvRll5H22KMYIiM9Xp5qaiL36mtoysvDNGIEw9auwWAyBXKVOl3j8ePkXrmE4pvu5GehEyioqGdgfAT/u+0vxDXVMfz993r0VZWm9RYislspNdWfeXjS2ipcRL4rIk87uln/p4j805+F9gTG2FjirlpC5dvvYK2sPG/85dtXE2Y18/SEpSCCzWDkmQuu5ZnxS6nZvJm8r30Nc1GRx8sr+++LNOXlkfj1W2g6fpzS558P5Op0iYo1a1BGIz+tSiO/oh4F5FfU80ziVJpOnqR227Zgh6hpWoB4chr4AjAAWAR8jL3rkOrODKq7SFi1CtXQQOW6tkUudXv3suDULtZmzCE/pm3v9LunLGDwM09jzjtJ3vUraDjacdsCS0kJJU8/TdTcOaQ88ACxV1xB6d+epTHX6xv5g0aZzVS++RZ708ZRaGx7xbU5ZSwVEbGUvfBCkKLTNC3QPEkeGUqpXwK1Sql/A1cC4zs3rO4hfPRoIiZOpPyll1vulFZWK4W/+Q3mxH6sG7eozfcjQo3ctyiL6LlzGfryyyDC6dvvwFxY6HY55556CltDAyk/vR+AlJ89gERGcvbBB3vMHdo1n3yCtaSEt1PPf46XxRDC20OnU7vlE5ry8ro+uG5EKUVPLirWtGaeJA+z42+FiIwD4oD0Touom0m4aZW9yMXRFLf81VdpPHSY9F8+wEMrppIWH4EAafER/P7a8SydZL8rPTxrJIOfew5bbS2nv/0dbLW1Tudf/+WXVLyxmsSvfY2w4cMACElOJuUn91G/azcVq4P/jBFls9Fw5IjbRFaxZg3G5GQKsiY5Hb93wiUQGkrZi3232a4ymzl18y0U/PhenUC0Hs+T5PGciCRg74b9LeAQ8IdOjaobiVm0CGNCAuUvvYylrIxzTz5F5MXTiVm8mKWT0th2/zxOPHIl2+6f15I4moVnjSTtySdoPHqU/B/fe14rLKUURQ//DmNCAsnf+XabcXHXXkvktGkUP/ZHLAF6bokvzEVFnL79Dk4sXcaZ7/8Aa/X5JZaWkhJqPvqYuGuu5p7LxxARamwzPiLUyJ3LphF7+WIq16zBWuM8kfZ2Jc8/T92uXVS99x5V774X7HA0zS8uk4eIpAAopf6ulCpXSm1RSg139Iz7rKvpehuDyUT88uXUbN7M2Qd+hq2ujgG/+AUizh5Bcr7o2bNJ+cXPqfnoI4r+0DbnVr33HvV79tjvgYiNbTNORBjwq/9BNTRQ+LvfBWx9vFH1/vvkXn0NdXv2EH/9cmo+/pgTy5fTcOxYm+9VvvU2WCzEX3stSyel8ftrxzu9Iku8+WZstbXn1SH1BQ1Hj1HyzN+IveJywi+YQNHDD2MpLw92WJrmu+Yy2PYvoBDYANwKxLn6XjBfU6ZMUV2h6cwZdWjUaHUoa5QqfPRRn+ZR+LvfqUNZo1Tpf/+rlFLKWlenjs29RB1ftkzZLBaX0517+ml1KGuUqtq0yafl+sJSWanO3HefOpQ1SuVev0I15OYqpZSq3bVLHZ01Sx2eOElVvP2OUkopm82mcq68Up1YcYNH885dsULlLL5c2azWTou/u7E1NancZdeqozNmKnNZmao/elQdGjdenbn3vmCHpvVRwC7l5/HXXa+6acB84Ebg9yKyHXgZeEspVd+5Ka17CU1LI2b+fOoPHCD52771Ft//Jz+h6dRpih7+HabBg6nftx9LYSFpf3wMMRpdTpd0221Uvfcehb/+DabBg7GUlGDOz3e8CjDn54MIkRddRNTF04mYMAEJDfV1Van9bAcFDzyApbiY5O99j+S77kRC7LtJ5JQpDFu9mvwf3UPBvfdSv38fsYsX05RznAG//pVH80/82s0U3HcfJ2/6Ggk33kDMokUYwsN9jrcrWKuqEJPJ5zhL//FPGg4dIv+HD/LNZ/dSUFHPXWMWcPXbbxN75RXEXHJJYAPWtC7g0U2CImICLseeSC4FNiqlburk2DrUmTcJtmerr0c1NmKMj/d9HrW15H3tZswnT6KsVmIum0fanzru5aVu715OrroJWv+vDAZCUlIITRuIqm+g4dAhUAqJjCTywqlETb+YqOkXETZiBNLBzYZNJ09Ss3UrtZ9speajjzANHcrAxx4lYsIEp99XZjPFf/wjZf/+DxIRAUqRufUTj7ofUTYb5S+8QPlLL9N08iSG2FjirrmGhBXXE5aZ2eH0XcVcWEj1hg+p3rCBul27MERFEX/tMhJWrsSUnu7xfBqzszlx7XVUTpnBNwYvpd5sr/cKtVr4y8dPkhpqZez693TXLVqXCsRNgh7fYS4imcBK4GvYm+06b1bThboyeQSKubCQvBU3YK2sZMT77xE6cKBH01Vv2oy1vJzQtDRCB6URmpLS5grDWlFB7eefU7t9O3XbP/uqSWxICKb0oYRlZBKWkWF/jRhO06lT9oSxdZu9+3QgdNAgYhcvIvk73/Ho7vjKd9/l7C9+SdySK53ehe+OUoq6z3dS8eqrVG/YgDKbiZg0ibhrrib60nmEpvTvcB7moiLq9+8nbPhwTMOHu62HstbUULP5I6o/+IDG3FxCU/oTkppKaOpAQlNTCR2YiiEmltrtn1K94UMa9u8HwJQxgpjL5mM+c4aq9evBYiFq9mwSblpF9Jw5bu+YVxYLeStXYT5zhu8v+ilHG9pe6I8qO8njW/5C0qobGfCg28fjaFpAdXryEJEhwA3Yk0YU8ArwilLqsD8LDZSemDwAzGfPYiktI2Lc2E5dRt2uXTQey6YxJ4fGnBzMZ860uXqRyEiipk0jatYsomfNJHToUI8bAjSzVlUh4eEed6Wybm8+j60/2tJ1yX2LslgyNILKdW9S8frrNJ2wd98ePmECMfPmEXPZPEwZGYgItsZG6nfvpuaTrdRu3UpjdnbLfI1JSUROnWp/TbuQsMxMbFVVVG/cRPUHH1D76acosxljv2QiJlyApeQcloKzTluyhY8bR8yCBcQsmM/7lWEt8Y42NfGA7SgDPn4Py7lzhA4eTMINK4hZtAjT4MHnzafk+ec59/ifSHviT1zwsQ1nv7Q7DrzJsuOfMPS/LxA51a/fss+UUjTl5GApLSMkOQljUhLG+Hiv9wXNO9aKCup276bu852IKZSEm28mtH/HJ02B0KnJQ0Q+xV7v8Tr2hNHtjtI9NXkEi62+nsbcXJpycwnp15/IyZM6LNIKpHV783lgzYGWohuwN+Ntbo2llKIxO5uaTZuo3riJhgMHAAgdOgTToMHU7d6NamhAQkOJmDqF6FmziJg0iabcXOp27qR2504sBWcBMMTG2u+tsVoJGZhK7IKFxCxaSMTEiW2uFmxNTViKijAXnMVaWkLEBRcQmpbmPt6rRzOv5BBl/32R+t27AQjLyiJm/nxiFswnLCuLpuPHObHsWqIvvZS0p55k1h82k19xflXhsCgDz3/0OBIayrA312EICwPsB3TL2bPUHzxIY04OprQ0wsdPwJQ+NCD9gymlaDx2jKr/+z+q/299S9JuERJCSGIixuQkjHFxYLGizGb7y2Kx/7VaMA0eQviYMYSPHUv42DGEpqW1STrWigoas7NpOHaMxmPZWMvKiLzwQqIvmYtpyBC/16MnsVZUULtzJ3Wf76Ru504ajx61FzWbTCirFQkJIeHGG0m6/VuEJCd3aiydnTzmAluUp+VaQaCTR9dwdrXQ/p4WT8x8ZJPTA2hafATb7p933nBzURE1mzdTvXETlsKzRF40nahZM4maNs1lsZo5P5/anTup370HY3w8MQsXEj5urE9n0Z7E23T6NNUfbqT6ww+p37MHlCJ00CAwGrBVVjH8nbcJSU52mzgX1J/k1K23EbdsGaGpA6g/eJCGg19iLSs7b9mGmBgixo8jfPwEIiaMJ3zUKEJSUz1KKLbGRhqzc6j+cIM9YeTlgcFA5IUXErt4EaZhw7CUlGItLcFSUoqltBRraan96jIkBAkNtReVhjrei9CYe4LGnBywWOzxxcYSPmYMEhJC47FjWIqLv4o9NhZjTIy9kQdgSk8neu5coufOIWLq1E7pCFQ1NdF4/DgNhw7TcPgwlpIS+7YyGMAgiHz1HpsCmw2lbGC1gbKhbAplMYPZYk+aLS8zoSkDiL1qCTFz57o8CVNKUb93L+UvvkTVBx+A2YyEhxMxaSKRF15I1LRphE+YgKWwkJJn/kblW28hoaEkrFpF0m23EpKUFPBtAl1c5+H1jO2dJy4BipVS4xzDHgOuApqA48A3lVIVjnEPALcBVuAHSqkOn8Skk0fn6+hqwRvD7n/XadGNACceudK/QDuBt/FaSkqo3ryZ6g8/pG7nLgb+7mFiFy9uGe8uCRf87OdUrlkDBgNhGRmEjxtH+LixRIwbR1hGBub8fOr3H6D+wH4a9h+w32vjOGBLRARhw4ZhGjGCsBHDMQ0bTkhiAk2nTtF43H6l2Zibay+2tNnsCWPaNGIXLyJm/ny/z3JtjY00Hsum4csvaTh0iIYvv0QpG+GZIwkbOZKwkZmEjRxJSP/+iIi9gcbHW6jZsoW6zz9HNTUhEREYo6PtvRjYbPa/VivYbIjJhDEhAWNiIiGJCRgTEu2f4+LsB/32rDaa8k7Q8OUhGrOzUWZ7JxmGyEhCBgz4KkE0JwvHMjEa2iQTMRhBxJ44Q0LsSTMktOVzw9GjWEtKMMbFEXvlFcRdcw3hEybYi1jr6qh85x3KX3qZxiNHMMTEELdsKbGLFxMxblxLsmm/T/xsQhSTPl5D5dvvIGFh9oYkI7MI6ZeMMSmJkOR+hCQm+NWiErp/8pgD1AD/aZU8FgKblFIWEfkDgFLqpyIyBnsz4GnAQOBDYKRSyu2DMQKdPAJ1ht2buDv7vm9Rllfby9srj0Dz9v/blfHamppoPHqMsIwRGCIiOv5+QwMNhw/TeCybptzjNB7PpTH3eEuxXTMxmTClp2MaMZywYcMxjRhO1PTpnXZG62obuxpuq6uj9rMd1G7fjmqoB4MRMRpADPaDucGIrbEBa1k51vJyrOVlWMrKsVZU2JOLC8a4OMLHjiF8zBjCRo8mfMwYTEMDU+TXTFks1H76KZXr3qR640ZUYyOmYcOImDiR6g8/xFZdTVhWFgmrVrFlyGQe/fhUm/UHXJ6YXR7XRMnTT1P13nv2xNZ+/eLjSfzGN0i+606fYu/WyQNARNKBd5qTR7txy4DlSqmbHFcdKKV+7xi3HvgfpZTbZ7sGMnkE8gy7N3F19g327ePN9grmNvZl2d1hn/A24dlqa2k8kYe1vAzT0KH2Ogg39xH5sgx383G2va6bksbq3fkB3Y7KZsNWV9e2+Xorhuhoj4oqA7Xu1upqqtevp3Ldm9Tv30/M/Pkk3LSKiMmTefOLAqfbJTzUQHmd+bx5tT45sdXXYyktxXLuHNbSUiwlJVjOlWApLSHq4hnELlrodawQmOQRzMfQ3gq86nifBnzWatwZxzBncd0B3AEwJIAVbo+tP9rmnwtQb7by2PqjfTp5DIyPcHr2bRTxens1Dw/G1Z0v/99gxgvnH4zzK+p5YM2BNrG1Z4iK8qoVny/LaJ6u/XZxtY1f3nEaa7uDvL+/LTEY/L43xtd1d8YYE0P88uXEL19+3jhX26X9sGYFrX5vhogITIMGYRo0yKt4uoK7O8ybHzWbBVyIvVNEsNdZbPFnoSLyc8ACvNg8yMnXnJ5SKKWeA54D+5WHP3G0VuDkAOlueKB11yKz+xZlOT1r8mTHd2bppLSgrJev/99gxQvuE17zeH/3F1+SqquDrqt9on3iaNZVvy1XOlr3QP0mvV3PgfER3fZ40JrL5KGU+hWAiHwATFZKVTs+/w/25rs+EZGvY69Iv6xVS64zQOuG8oOAAl+X4QtXZ9gD4zsuf/ZXIM+A/InB2c7q6uz7sfVHg7a9fBHo/29X/LhdHXTaH6z92V98SaquDrpGEaeJwtXwYO8r7tY9kL9JV/tefEQojRbbeSdml47qF/TjgSc8qT0agr11VLMmfHyeh4gsBn4KXK2Uqms16i3gRhEJE5FhQCbwuS/L8NV9i7KcdiXeXLHVmTo6w/TWur35zHxkE8Puf5eZj2xi3d78Dr//wJoDbR4d+8CaAy3TOet6PtDby9uYvZ2mo3i9mVdH2ytQXB1c3RUZBmoZ7g7srg66VqWcbuOVFw32aV/xZZ/whrt1D+Rv0tW+9z9Xj3XaA/XmI+fcLruzt4unOqwwdxQxrQDWYi9KWga8ppRy20+4iLwMXAIkA0XAQ8ADQBhQ6vjaZ0qpu1ot51bsxVk/VEq931HwPleYv38/FB44b3BJTSOnyupostowGQ0MSYwkOTrM+/l7Oa/PTpQ6HQ4wfZh3rWJKahrJLanF1ur/ahBheHIUydFhTuNq/tyeyWhg8pAEt8sKxPbqKGZvpwGcxuUqXm+Xv+dUuU/by1uu4rK5+c0Gcn8B59vR3fq33p882fa+xOXrb9KbZeScq3E5XUa/aK/3+0AdDzL6RbfEfMg2lF9bbvGp8UGnVpg3U0o9LCLvA7Mdg76plNrrwXQrnQz+h7vlAA93NN/OlBwdFpAds/1O2WS1kVtS27KM9kxGg8sfo7dOldWdd3CxKcWpMvuFnrO4XB2MnMXUWqC2l7uYXc3f1TR5pbX25vsutr2z+Xm7fFfbpaPt5a3mZXuT7AO1DHC+rwAMSYx0etBtPiA622bu9hVXJzTe7hOu5tVR8nI23NU2DjGIV79tT9a/PXfHA2fbJVgNezpMHg6RQJVS6v+JSD8RGaaUOtHhVN3V5Y906uyveWQT+Y1O7g+ojWDb98+/P+CUy24wxjPZyx3iOjc3tg2sjXAal6sy6bT4CLZ9s/Pvv3AX8xMTJzqtW3A1jSuutn1Hyz/xzfNvBvyxu/s/Ary9kh2v1rYGcH9xtYyZ7vbh++ex1Umdz7RANe8tdt0gw90+4Wpe1w1OY3Vu/nnDfz/Lfrbu7BZJV7/JcIOL5rVu9i9vuTse/OjVL5zuq8FofNBh8hCRh4Cp2Ftd/T8gFPgvMLNzQ+u5vK2EDGSTUHcVwx2VVbffWbuivqc5Nmcxx0WEuqw4dDWNK+5+XN5WprtqgdZV28vX/cWbSv6O9uFAtULztvLd3T4RyKbCrrbxj179wun3myvZA/Ebdvf/7U4NVTy58lgGTAL2ACilCkQkxv0kfZsvLXt8+TE621ndHdhc7Xi+3C0eSK5iFsFlxaGraVzdeOVu23ubDIJ9/0dzDN4sz9vWQ13V+tDbExp3+4S7eXmz7GbOtrGr35C7pOZrAnE2XbBPXFrzpJC0ydGkVgGISFTnhtS5uqKlQle03HLV4gdw+Qxxd3E5a1HVVVw997zCSRIA+4/e1TQPXTXW623v7rnr7qYJ1Pbqin3S29ZDXdX60FUyav4feLNPuGud5s2y3XG1XexJzYIh/AyhCduRkAq/Wky60ryvDow3IoZGj/bVzuLJlcdrIvIsEC8it2NvEfX3zg2rc3TV/RSBPjP15m7ex9YfdXkw6w5nzK54c5bX/KN3d/bt7ToG8mZAb4ovArFPmq1mTladJLcyt+VV0VDBkNghDIsbxrC4YZytKQRiaX8/blcUpbrTciZtaQCxgjIQERrCPQtHcPXEVJZOSqPR2kh5QzllDWX0759HaUM5YqxFjPU030scEx7K5BFJVBwpxmxV2MyJ2BrSCLMN4Lopg512j+JLInS2Xe64NJGHt7xI5IA9GMPsvQirlLew1IymqOxilLrUp16dz9ac5d0T71JcV0xJfQml9aWUNpRSUl9CbWotP1p4Oz+YfK3X8w0UTx9DuwBYiH3PW6+U2tDZgXnC26a6we6Yzxeu+gtyV6HYHXuo9UV36FvKW+v25vPA2xuxmHIwV10AtnC3MTfvk4bwU5gSt2FrHIC1YRApphF8+tOrzvt+rbmWQ6WHOFBygIMlB8kuz+Z09WmsrfoQHRg1kITwBE5VnaLaXN0yXFlN2Jr6Y6nNxFI1HltjKmnxkQHd92vNtRwpO4JN2ciIzyAh3HnTZaUUuZW5bDmzhbVHPuREzUEQ71urKWU/KLfuXFe1qlIOkVCyEkcSoYZw8EQ05ZVRJEaHsmzSQKYNT0AphVVZUSgMGDCIARHBgAGjwYggGA1GQg2hhBhCCJEQQo2hhEgIh8sO89bxt9hxdgcKhaVuKJbKyVjrhxIS+wWh8TsxhNSSHpvOjaNu5OoRVxNj6rjEv7C2kL8f+Durs1djsVmIMcWQHJFMckQySeFJ9r8RSUzuP5nJKZO93mbQRR0jisgflFI/7WhYMHibPHpal+DgOuG5bSHVTROhL9buOcMfNm2hpL6Y/qYsfrJwgkeJo6apBpPRhMnY+Q+7UkpxsOQgG09t5P998Ta2EPvZp7lyIg0FNwKu/y/D7n8XJU1EDX8SCalCDJaWcYOiBzE2eSwZ8Rnk1+RzsOQgxyuOtxwcB0UPYlTiKIbFDWN4/HCGxw3nYF4YT204SUFFPanx4Xx7Xn9GDWlg7cE9vPnlF9hC8zFG5iGiUE3JzB00n+9fdB1ZCVltzo6VUlQ1VVHWUEZ1U3XLAbT5IBpqCEUQTlSd4HDpYfur7DB5VXlt1i85IpnM+EwyEjLIjM8k1hTL9rPO6kXLAAAgAElEQVTb+eTMJxTU2juRyEzIZFbaLBLDElsO5FabFZuyYcOGyWAiITyBhPAEEsMT2Z3bxD8/LuZsOQyMjzrvishqs3Ky+mRLXEfKjnCo7BDVTdUEWlp0GleNuIrwhmn88d3Stic6Jhs3XFJGdv0H7C/ZT0RIBDMGzuDCARcyNWUqmQmZGOSrmoPiumL+ceAfvH7sdRSKZRnLuH387aRGpwY87q5KHnuUUpPbDduvlJrgz4IDoS9ceQSyV9uewmqzsu/cPjad2sTm05s5VX0KgFBDKJP7T2ZG2gxmDpzJyISRLQe80vpS9hTvYVfhLnYV7SK7PJvI0EjmpM1h/tD5zEqbRWRox89l98bBkoO8mfMmm05voriuGKMYaaxJx1I9FkNoGaakrdSfvhlLzViXJygzH9lEiel1TElbqcu7A2vjAIzhBSQkFDFnfAOHSg+RX5NPQlgC45LHMT55POOSxzEuedx5Z/UdXam1FKdVnyO5fzZpg45xonYfNmVjSMwQhsYOpbShlNL6UsoayjDbnNcvODMwaiCjk0YzOnE0o5NGEyIhZFdkk12eTXZFNscrjtNobbTHFBLB9NTpzB40m9lpsxkQNcC3f4AXlFLk1+RT1lCGUYz2qwsxIHz1V6GwKZs9eSlry1WJ1WbFYrPYX8qC2WbGYrPQP7I/F/S7oCUBuCuu/LL0S9449gbbC7aTX2Ov04oLi2Ny/8lcOOBCCmoKeP3Y61hsFpZmLOX2CbeTFt15v+POfpLgt4HvACOAnFajYoBPlVI3+bPgQPA2efTEYhB/n6dhtVnJr8knpyKH4xXHyanI4UTlCfpH9mfGwBnMTJvJkJghQX9etdVmZVvBNjac3MCWM1soaygjxBDCRQMu4tLBl5IWk8aOszvYVrCN7HL7s8uTI5KZ2G9iSzk/2A9MF/S7gMn9J1NUV8SmU5sobywn3BjOzLSZLBi6gAn9JlBaX8rZ2rMU1BRwtvYshbWFlDWUcUG/C5g/dD4T+03EaDi/K/MmaxPr89bz8pGXOVBygDBjGDMGzuCyIZcxd9Bcrnxyj+P/ZSFy2F+RkBpqj/+ItNgkpyco/7ttI89m/whzxTQaC5fZ16HdPllnriMiJKLD/5EvJ0dlDWVsPLWRD09+SFlDGUnhSSRFJLX8TQxPJC4sDqvNitlmbjlwmm1mrDYrg2MHMzpxtMviqWZWm5XT1acpbyxnbNLYLrki7K4KagrYVbSLXYW72Fm4kzM1ZzCKkatGXMUdE+5gcMzgjmfip85OHnFAAvB74P5Wo6qVUuc/HzMIfOmepDv0Vvmvg/8iryqPhy5+qMMDgrcJr9Zcy+6i3Xx29jN2F+1uc8YHkBqVyrC4YZyuPs3p6tOA/dJ75sCZzEibwfTU6USFetagzmw188LhF5jUfxKT+k/yaJr2CmsLWZu9ljU5ayisLSQmNIZZg2Yxb8g8Zg2cRbTp/G63i2qL+LTgUz4t+JQDJQcYFjeMqSlTmZIyhbFJYwk1fvWUNYvNwp6iPWw4uYGNpzZyrv7cefOLNcWSGpVKjCmG/ef202RrIjE8kXlD5jF/yHymDZhGaUMprx19jdXZqylrKCM9Np2Vo1Zy9Yir28TY+v9lCMsncthfUdWTeHjWw+f9v8xWMyveWUFxbQXq9H2cLVd+7ZM9sVhWs/8GBCElKqXLltlVxVbTgS9b9aobA4xRSu3wZ8GB0BMfQ7s2ey0PfvogAI/NfYzF6Ys7mALe2H2C3+7+MU2qijCSmJI2glnpIxkYPZC06DTqLfXsKNzBjrM7OHDuABZlwWQwMbH/REYljiIjPoMR8SMYET+iTWI4XXWabQXb2Fawjc/Pfk6dpY5YUyy/nP5LFg9zH1dxXTE//ujHfHHuC0Ynjua1q17zeBtYbVa25m/ljWNvsCV/C0opZgycwfKRy5k7aG6bg38g2ZSN/ef2c7ziOClRKaRGpTIgakCbbVJrruWT/E/YeHIjW85soc5SR3RoNPWWehSKuYPmsnLUSqanTneZ+FufoCQO2kxTzHr+Mu8vzB08t833ntn3DE9/8TR/veyvzBk0x+/164nFslpwdFXy2Iu9S/bm+zwMwK729SDB0JXJY1v+NvJr8lmRtcLneXx+9nNu/+AOpDGDJlVFSEgDv5j4/1gxZYTb6Z7+4mme2fcMM9NmUlpfSn5N/nmVfwYxMDZpLBelXsT01OlM7D+RMKPn/U6ZrWa+OPcFT+5+kv0l+7ly+JX87KKfEWuKPe+7Owt3ct/H91FnqWPmwJl8eOpD3rjqDbISO276WFhbyDf+7xvk1+STHJHMsoxlXJt5LYNiut/DbhqtjXxW8BmbT28mPiye67Ou97oc2mw1c8O7N1DZUMmaa9YQFxYHQE55Dte/cz0Lhi7g0TmPBiTenlgsqwVHVyWPL5RSE9sN65EV5r6y2qxcufZK8mvy+e8V/+WCfhd4PY/cylxuePsm6uuiqMn7NsbwAiKHPoe1dCG/u/Qelz/u3Ipcrnv7OhYOXcgf5vyhZXh1UzUFNQUU1NhbrExOmdxyYPKHxWbh+f3P8+z+Z+kX2Y+HZz7MtNRpgL3S8T+H/sMTu59gcMxgnrjkCZIjkpn3+jxuyLqBn07ruAHe47se54VDL/DonEe5dMilhBo65yqjO/my9Etuevcmlgxfwm9n/Rarzcot79/CqepTvLn0TRLDEwO2rO5QLKt1f13Sqy6QKyI/AJ5xfP4OkOvPQnuarflbya/JJ8QQwsOfPczLV77stDLVlfKGcr774XdpbILa018HWzjWuuGYq8YTkrCZP2yYwdJJ1503nU3Z+NX2XxEVGsVPLvxJm3ExphiyErM8Otv3RoghhG9P/Daz0mbxwNYHuO2D27hlzC18a/y3eHjHw6zPW89lQy7jtzN/21LWP2/IPN7JfYcfTfmR24rQOnMdq7NXM3/ofBam+/bs5Z5obNJYbh13K88feJ5F6Ys4WXWS/SX7eWT2IwFNHBDcJx9qfYsn3ZPcBcwA8rE/8e8iHM8Q7yteOfoKyRHJ/HrGrzlcdpjV2as9nrbR2sjdm+/mXP05ak/dgjJ/dbBoLL4cUJSHr3U67ers1ewp3sO9U+8lKcK75zT4a3y/8by25DVuyLqB/xz6D/Nen8eGkxv44eQf8sQlT7SpJF6asZSKxgo+Ov2R23m+e+JdqpuqWTVqVSdH3/3cdcFdjIgbwUOfPsSf9/6Z2WmzuWLYFcEOS9N81mHyUEoVK6VuVEr1V0qlKKVWKaWKuyK47uB01Wm25W9j+cjlLBm+hGkDpvHUnqcobyjvcFqlFA9ue5C9xXt5eNbDDAhve5WgzIk0lc4hNG4fe4r2tBl3ru4cT+x6gmkDpnHNiGsCuk6eigyN5BfTf8HTlz3NlP5TeHbBs9w2/rbzKoovTr2Y/pH9WZezzuW8lFK8dPglRieO9rllVk9mMpr47azfUtpQiiA8ePGDQW8erWn+6DB5iMhIEdkoIgcdnyeIyC86P7Tu4bVjr2EQA8szlyMiPDDtAerMdTy156kOp31639O8d+I97p58N4vSFzntVM1YdRmxock88vkj2NRX3TM88vkjNFobu8VBZvag2fx90d+Znjrd6Xijwcg1I65hW8E2imqLnH5nZ+FOcipyWDV6VdDXJ1jGJY/j0TmP8ud5f+6SG+M0rTN5Umz1PPbHx5oBlFL7gRs7M6juosHSwNqctWTFXsy1//slw+5/l68/e4ppSVezJnsNB0sOOp3OYrPw+x2/52/7/sbSjKXcNu42wEXvrcum8PPp93G47DBv5rwJwEenP+KDkx9w1wV3MTR2aFetrl+WZizFpmy8nfu20/EvHn6RhLAELh92eRdH1r0sSl/ERakXBTsMTfObJxXmkUqpz9udLVpcfbmZiPwTWAIUK6XGOYYlAq8C6UAesEIpVS72mT8FXAHUAd9QSu1xNt+utD5vPZWNlZTljqWuyt5+Pr+inrIdk4jP2szDnz3Mi1e+2KZ/mpqmGu7bch9b87dyy5hbuGfKPW3OtJ1VaCo1kFeOvsKTe55kxsAZPLzjYTLiM/jG2G90yXoGwpDYIUxJmcLa7LXcNq5t0VZ+TT4fnfmI28bd5lXzYU3Tui9PrjxKRGQEXz3PYzlw1oPp/gW0v9PsfmCjUioT2MhXd65fDmQ6XnfwVcuuoHrlyCsYzCnUVaW3GV7fFIr13BIOlh5kTfaaluEFNQXc/P7NbC/Yzi+n/5L7LrzPo1ZZIsJPp/2U8oZyVr67kqLaIh66+KFOu1musyzLWMap6lPsLW77iPtXj7yKIH7dI6NpWvfiSfL4LvAsMEpE8oEfYm+B5ZZSagvQvhuTa4B/O97/G1jaavh/lN1n2J8dEviuJL1wsOQgB0sPUl96Ee2fgQBwrnAMU1Km8NSep6hsrGT/uf2sencVRbVFPD3/aa8PlGOTxrIscxnn6s9xQ9YNTOw/seOJupkFQxcQGRLJ2pyvWo/VW+pZnb2ay4Zcpsv5Na0X8aS1Va5Saj7QDxillJqllDrp4/JSlFJnHfM9C/R3DE8DTrf63hnHsKB59eirRIREkMwMp+MHxkfys4t+RnVTNXdvvptb199KeEg4L1zxAjMGOp+mI/dMuYe7J9/ND6f80J/QgyYyNJLFwxazPm89deY6AN7NfZeqpipWje57zXM1rTfzpLVVkoj8GfgE+EhEnhKRQN904Kz5jdNb30XkDhHZJSK7zp07v5O7QKhoqOD9E++zZPgSfrJwosvHcY5MGMnKUSvZXbSbMUljeOnKlxgR776rEXfiwuL41vhvedwxYXe0LGMZ9ZZ61uettzfPPfISoxJHMbl/0Huz0TQtgDypMH8F2AI03wJ9E/ZK7/k+LK9IRFKVUmcdxVLN94ucAVr3QzwIKHA2A6XUc8BzYO+exIcYOvTm8TdptDZyQ9YNZCW6fxznDyb/gLHJY1kwdIGuDAYu6HcB6bHprMtZx6CYQWSXZ/PrGb/us81zNa238iR5JCqlftPq829FZKnLb7v3FvB14BHH3zdbDf+eiLyC/Q72yubira5mUzZePfoqk/tPbun6w12XDxEhESwZvqQrQ+zWRISlGUt5cs+TPL7rceLD4vt881xN6408qTDfLCI3iojB8VoBvNvRRCLyMrAdyBKRMyJyG/aksUBEsoEFjs8A72HvLysH+30l3/FhXQLi04JPOV19mhtH9YlbWTrF1SOuxihGviz9kusyryM8JDzYIWmaFmCeXHncCdwDvOD4bARqReQeQCmlzu+z2z5ipYv5Xebkuwp7q66ge/XIqySGJzJ/iC+lchpAv8h+zEqbxdb8rdyQdUOww9E0rRN0mDyUUjFdEUh3UFhbyJb8Ldw27rYed49Fd3P/tPs5UXmC1OigtrjWNK2TeNLa6rZ2n40i8lDnhRQ863LWYVM2rs28Ntih9HiDYgYxe9DsYIehaVon8aTO4zIReU9EUkVkPPAZ0OuuRmzKxrqcdVyUelG3fKqdpmlad+JJsdUqEbkBOIC936mVSqltnR5ZF9txdgf5Nfn8YNIPgh2Kpmlat+dJsVUmcDewGntnhjeLSGQnx9Xl1mavJdYUy2VDz6vP1zRN09rxpNjqbeCXSqk7gblANrCzU6PqYhUNFXx46kOWDF+ib/TTNE3zgCdNdacppaqgpUnt4yLyVueG1bXePfEuZptZV5RrmqZ5yOWVh4j8BEApVSUi17cb/c1OjaoLKaVYk72GsUljW+4o1zRN09xzV2zV+hbrB9qNa/+cjh7rUOkhjpUf01cdmqZpXnCXPMTFe2efe6zV2asJN4br/pc0TdO84C55KBfvnX3ukeot9bx/4n0Wpi8kxtTrbl3RNE3rNO4qzC8QkSrsVxkRjvc4PveKnu42nNxAjbmGZRnLgh2Kpmlaj+IyeSilOn74dg+3+thqhsYOZUrKlGCHomma1qN4cp9Hr5RXmcee4j0szViqH1SkaZrmpT6bPNbmrMUoRq4ZcU2wQ9E0Tetx+mTyMNvMvJnzJrMHzaZfZL9gh6Npmtbj9Mnk8cmZTyhtKOXaDH1vh6Zpmi/6ZPIYmTCSOyfcqZ83oWma5iNP+rbqdQbFDOJ7k74X7DA0TdN6rD555aFpmqb5RycPTdM0zWti72W9ZxKRc8DJDr6WDJR0QTjdUV9ed+jb69+X1x369vp7su5DlVJ+NTXt0cnDEyKySyk1NdhxBENfXnfo2+vfl9cd+vb6d9W662IrTdM0zWs6eWiapmle6wvJ47lgBxBEfXndoW+vf19ed+jb698l697r6zw0TdO0wOsLVx6apmlagOnkoWmapnmtVycPEVksIkdFJEdE7g92PN4QkX+KSLGIHGw1LFFENohItuNvgmO4iMifHeu5X0Qmt5rm647vZ4vI11sNnyIiBxzT/FkcDzVxtYyuJCKDRWSziBwWkS9F5O6+sv4iEi4in4vIPse6/8oxfJiI7HDE9aqImBzDwxyfcxzj01vN6wHH8KMisqjVcKe/C1fL6GoiYhSRvSLyjru4eum65zn2yy9EZJdjWPfc75VSvfIFGIHjwHDABOwDxgQ7Li/inwNMBg62GvYocL/j/f3AHxzvrwDex/6I4OnADsfwRCDX8TfB8T7BMe5z4GLHNO8Dl7tbRheveyow2fE+BjgGjOkL6++IJ9rxPhTY4Vin14AbHcP/Bnzb8f47wN8c728EXnW8H+PY58OAYY7fgtHd78LVMoLw/78HeAl4x11cvXTd84DkdsO65X7f5RunC/8JFwPrW31+AHgg2HF5uQ7ptE0eR4FUx/tU4Kjj/bPAyvbfA1YCz7Ya/qxjWCpwpNXwlu+5WkaQt8ObwIK+tv5AJLAHuAj7HcMhjuEt+zawHrjY8T7E8T1pv783f8/V78IxjdNldPE6DwI2AvOAd9zF1dvW3bHsPM5PHt1yv/ep2EpEDCKywpdpu1AacLrV5zOOYT1ZilLqLIDjb3/HcFfr6m74GSfD3S0jKBxFEZOwn4H3ifV3FNt8ARQDG7CfLVcopSxO4m1ZR8f4SiAJ77dJkptldKUngZ8ANsdnd3H1tnUHUMAHIrJbRO5wDOuW+71PyUMpZQO6e5/mzh5M3lvbJbtaV2+HdysiEg2sBn6olKpy91Unw3rs+iulrEqpidjPwqcBo519zfE3UOse9G0iIkuAYqXU7taDnXy11617KzOVUpOBy4HvisgcN98N6nr6U2G+QUTuFXvlZmLzy4/5BdoZYHCrz4OAgiDFEihFIpIK4Phb7Bjual3dDR/kZLi7ZXQpEQnFnjheVEqt6SC2Xrf+AEqpCuAj7OXZ8SLS/Pyd1vG2rKNjfBxQhvfbpMTNMrrKTOBqEckDXsFedPWkm7h607oDoJQqcPwtBtZiP3nolvu9P8njVuC7wBZgt+O1y4/5BdpOINPRisKEvULtrSDH5K+3gOaWE1/HXhfQPPwWR+uL6UCl49JzPbBQRBIcrScWYi/LPQtUi8h0R2uLW9rNy9kyuowjpn8Ah5VSf2o1qtevv4j0E5F4x/sIYD5wGNgMLHcSV+t4lwOblL3g+i3gRkeLpGFAJvbKUqe/C8c0rpbRJZRSDyilBiml0h1xbVJK3eQmrl6z7gAiEiUiMc3vse+vB+mu+30wKoW6sPLpCuwtdY4DPw92PF7G/jJwFjBjP2O4DXvZ7EYg2/E30fFdAf7qWM8DwNRW87kVyHG8vtlq+FTHjnkc+Atf9TbgdBldvO6zsF9O7we+cLyu6AvrD0wA9jrW/SDwoGP4cOwHwBzgdSDMMTzc8TnHMX54q3n93LF+R3G0qnH3u3C1jCDt/5fwVWurPrHujhj2OV5fNsfXXfd7n7sncRQrfBt7k1KwX14/q5Qy+zRDTdM0rcfwJ3n8HXs79H87Bt0MWJVS3wpQbJqmaVo35U/y2KeUuqCjYZqmaVrv40+FuVVERjR/EJHhgNX/kDRN07TuLqTjr7h0H7BZRHKxV9wMBb4ZkKg8lJycrNLT07tykZqmaT3e7t27S5SfzzD3KXmIiAGox94ELgt78jiilGr0Jxhvpaens2tXd2odrGma1v2JyEl/5+FT8lBK2UTkcaXUxdibFGqapml9iD91Hh+IyHXNXfr2JCU1jfzfwbM0WnQVjaZpmi/8SR73YL+ZplFEqkSkWkTc9T/UbXyWW8pd/93D8eLaYIeiaZrWI/naq64AY5VSBqWUSSkVq5SKUUrFBji+TpGVEgPAsaLqIEeiaZrWM/naq67C3mlXj5SeHEWoUTiqk4emaZpP/Cm2+kxELgxYJF0o1GhgeHI02Tp5aJqm+cSf+zwuBe50NPmqxd5cVymlJgQksk6WmRLNvjMVwQ5D0zStR/IneVwesCiCICslhnf2n6WuyUKkyZ/NoGma1vd4XWwlIvMAlFInAYNS6mTzC5gS6AA7y8gB9krz7KKaIEeiaZrW8/hS5/HHVu9Xtxv3C09mICL/FJFiETnYatj/iEi+iHzheF3hQ2weG+locaUrzTVN07znS/IQF++dfXblX8BiJ8OfUEpNdLze8yE2jw1JjCQsxKArzTVN03zgS/JQLt47++x8Bkptwf6s4aAxGoSM/tEc1cVWmqZpXvOlpni4iLyF/Sqj+T2Oz8P8jOd7InIL9meh/1gpVe7n/NzKSolhe25pZy5C0zStV/IleVzT6v0f241r/9kbzwC/wX718hvgcezP4W1DRO4A7gAYMmSIH4uzV5qv2ZtPZb2ZuIhQv+alaZrWl3idPJRSH3dGIEqpoub3IvI88I6L7z0HPAcwdepU3x6D6DAyJRqA7KJqpqYn+jMrTdO0PsWfO8wDSkRSW31cBhx09d1AGdnSx5Wu99A0TfNGUO6OE5GXgUuAZBE5AzwEXCIiE7EXW+UBd3Z2HGnxEUSZjLqDRE3TNC/5nTxEJEop5VXf5kqplU4G/8PfWLwlImSmxHC0UCcPTdM0b/hcbCUiM0TkEHDY8fkCEXk6YJF1kayUGLKLdfLQNE3zhj91Hk8Ai4BSAKXUPmBOIILqSpkp0ZTUNFFa06WPX9c0TevR/KowV0qdbjeoxz3XNWuArjTXNE3zlj/J47SIzACUiJhE5F4cRVg9yUj9VEFN0zSv+ZM87gK+C6QBZ4CJjs89Sv+YMOIiQnUHiZqmaV7wqbWViBiBm5VSNwU4ni4nIvZKc508NE3TPObrM8yttO2mpEfLTInmaGE19keza5qmaR3xp9hqm4j8RURmi8jk5lfAIutCWQNiqGqwUFytW1xpmqZ5wp+bBGc4/v661TAFzPNjnkGR2d/xYKjCalJiw4McjaZpWvfnc/JQSl0ayECCqbmDxGNF1cwZ2S/I0WiapnV/fnVPIiJXAmOBltN1pdSvXU/RPSVFh5EcHaab62qapnnIn+5J/gbcAHwf+4OgrgeGBiiuLjcyRT9VUNM0zVP+VJjPUErdApQrpX4FXAwMDkxYXW9kSgw5RdXYbLrFlaZpWkf8SR71jr91IjIQMOP/Y2iDZmRKDLVNVvIr6jv+sqZpWh/nT/J4R0TigceAPdifwfFKIIIKhqwBX1Waa5qmae75nDyUUr9RSlUopVZjr+sYpZT6ZeBC61qZ+qmCmqZpHvO5tZWI3OJkGEqp//gXUnDEhoeSGheurzw0TdM84E9T3QtbvQ8HLsNefNUjkwfY6z30UwU1TdM65s9Ngt9v/VlE4oAX/I4oiEamRLM9txSrTWE0SLDD0TRN67b8ehhUO3VAZgDn1+VGpsTQZLFxstSrR7Jrmqb1Of7UebyNvS8rsCehMcBrgQgqWFo/VXB4v+ggR6NpmtZ9+VPn8cdW7y3ASaXUGT/jCaqM/tGI2DtIXDxuQLDD0TRN67b8qfP42NdpReSfwBKgWCk1zjEsEXgVSMd+z8gKpVS5r8vwRaQphKGJkRwprOrKxfZKh89WceBMJSsu7LGdDmia5oY/fVtVi0iVk1e1iHR09P0XsLjdsPuBjUqpTGCj43OXG50ayxHd4spvz3+Sy/1r9tNgtgY7FE3TOoE/FeZPYD/ApwGDgJ8Cv1VKxSilYt1NqJTaApS1G3wN8G/H+38DS/2IzWejBsSSV1pLXZMlGIvvNbKLarApdNNnTeul/Ekei5RSTyulqpVSVUqpZ4Dr/JhfilLqLIDjb39nXxKRO0Rkl4jsOnfunB+Lc250agxKoa8+/GCzKXKK7XfqHzqriwA1rTfyJ3lYReQmETGKiEFEbgI6vYxCKfWcUmqqUmpqv36Bf3DT6FT7RdORszp5+OpMeT31juKqQwU6eWhab+RP8lgFrACKgGLsz/NY5cf8ikQkFcDxt9iPeflsUEIEMWEhHNZnzD5r7uIlOixEX3loWi/lT8eIeUqpa5RSyY7XUqVUnh+xvAV83fH+68CbfszLZyLCqNQYnTz8cKzYnjwWjR3A4bNV+hkpmtYLeZ08ROR2Ecl0vBcR+aeIVIrIfhGZ7OE8Xga2A1kickZEbgMeARaISDawwPE5KJpbXCmlD3q+yC6qITUunIuGJ1LXZCVP37Gvab2OL/d53I29qS3ASuACYDgwCXgKmN3RDJRSK12MusyHeAJu1IBYahpPcqa8nsGJkcEOp8c5VlRNZkoMYxz1R4fOVuk79jWtl/Gl2MqilDI73i8B/qOUKlVKfQhEBS604Bmdau+mxJPy+gazldW7z2Cx2jo7rB7B6mhpNbJ/NJkp0YQYRFeaa1ov5EvysIlIqog0d8P+YatxEYEJK7iyBsQggkf1Hmv35vPj1/fx2q4e3TNLwJwuq6PRYiMzJZqwECMZ/aN1pbmm9UK+JI8HgV3YuxB5Syn1JYCIzAVyAxda8ESaQhiWFOVRc92tOSUA/O+mbH03NV+1tGp+MuOYgbH6ykPTeiGvk4dS6h3sj50drZS6vdWoXcANgQos2EalxnC4gz6ubDbF9uOljOgXxdnKBl7acaqLouu+sh03ByZ1p58AAB8HSURBVGb2t9dxjEmNpbi6kXPVjcEMS9O0APOpqa5SytK+00KlVK1Sqtc8AHz0gFhOltZR0+i6m5IjhdWU1Tbx7UsymDEiiac/yvGrWxOlFLnnevYmzC6qZmBcODHhoYD9ygM8KwLUNK3nCOTDoHqV5jvN3fXN9Olxe5HVzIwkfrxwJCU1Tfzr0zyfl/nRsXPMe/xjdp/s0s6EA+pYUU1LkRXQpsWVpmm9h04eLoxytLhyd8a8NaeE4clRpMZFMGVoIpdm9ePZj3OpajC7nMadj4/a++raeLjIp+mDzWpTHD9Xw8iUr5rlxkeaSIuP0PUemtbL+JU8RCRNRGaIyJzmV6ACC7a0+Ahiw113U9JksfH5iTJmZiS3DPvxwiwq6838/ZMTPi3zs9xSALZkB77Dx65wqqWlVUyb4aNTY/WVh6b1Mv48hvYP2CvID/FVh4gK2BKAuILO3k2J62d77DtTQV2TlZkZSS3DxqXFcfm4Afxz6wm+MSOdxCiTx8srrWnkSGE1ydEmDuZXUVLTSHJ0mN/r0ZWaW1qNbJc8xgyMZdORIuqbrESYjMEITdO0APPnymMpkKWUukIpdZXjdXWgAusORg+I4YiLvpm25ZQgAtOHJ7UZ/qMFI6ltsvDsx8e9WtaOE/bHm3x/XiYAW7NLfIw6eLKbm+n2b3s3+ZjUWPuzPYp0T8Wa1lv4kzxygdBABdIdjU6NpbbJyunyuvPGfZpTyriBccRHtr26GJkSw9KJafx7ex7FVQ0eL2v78VKiTEZunDaYxCgTW471vKKrY0U1pMVHEBXW9oJ2rKPFla730LTew5/kUQd8ISLPisifm1+BCqw7aG5xdbjdzYK1jRb2nCpnRkaSs8m4+7JMzFbFXzfneLys7bmlXDgskbAQI7MyktmSXdLjeqM9VlTdprK8WXM394fOVgYhKk3TOoM/yeMt4DfAp8DuVq9eY2RKDAYn3ZR8nleGxaaY1aqyvLX05ChWTB3ES5+f4oyTq5b2iqsayCmu4WJHEdickf0oqWnsUZXMFquN3HO151WWg73+aLS+01zTehV/nufxb2evQAYXbBEmI+nJUeclj09zSjAZDUwdmuhy2u/Ny0QQj1pebXe0srp4hCN5ZNqTUk9qdXWyrI4mq+28+o5mYxyND6w97GpK0zTnfE4eIpIpIm+IyCERyW1+BTK47mC0kxZX23JKmTw03m3LobT4CBb8//bOPD6q6l7g35PJRiAbJIGEJGwhCWEJgaAsKosL7kVRFOtSl/p86qO2VSvV0mpr7ROea+0TfNaqFeouWpVNIcoiyBYgkIQkJCRkD2QnySzn/XHvhARmkkxmkhmS8/185nPvPXf7/WbO3N895/c7v5M0lLX7T9Bi6jjj7g95VQT6ezM+KhiAiCB/EocFnld+j6N2Iq2sjI8KUnN7KBR9CGe6rd4C/hcwAXOBd4B3XSGUJ5EUGcTxk43U6QP/Tja0cLikllljbHdZteWmqdGcajTybWbHM+ruyK3iwlGDMXiJ1rLZ8eHsKThFQwfpUTyJ7DItrUqcvZaHcporFH0KZ4zHACnlN4CQUhZIKf8AzHONWJ5D4jDtTdqapmRHrtbFNGts58bj4rFhhAf68dEe++naS2pOk1/VeE7I7+z4cIxm2Xo/Tye7rI7o0HMjrayMjQjExyDOKz9OT7B8fSZvbeveIFKFwpNwxng0CSG8gKNCiIeFEDcAES6Sy2NojbjSjcfWnEoC/byZNDy403O9DV7cmDKczVnldrPKWo2D1d9hZerIUAb4GDzG79FZ5FdOeb3dLisAX28v4iIC+3XLo7HFxBvfHWNlWp6a4lhx3uOM8XgECACWAFOB24G7XCGUJxEZ7E/wAJ9Wp/n23EouHD0Yb0PXvrqFU6MxWyRr95+wuX9HbhUhAT6MGxbUrtzP28CMMUM8wu+xPaeS1Gc3tSaCPJszkVYdTzWb1M/TlGzLqaLFbKG0tqm1m0+hOF9xJtrqRz0F+ykp5d1SyoVSyh9cKJtHIIQgcVggR0pqKTrVSEFVIzO74O+wEj80kOToYD7aU2TzbXNHnubv8Grj77Byydgw8qsaKXCjk7nJaGbppwc52dDCM18cthktlV+lRVrFR9hveYDm96ioa6a8ruuDJ/sS32aW4+et/eW2ZHXsB1MoPB1noq1mCCEOA0f07WQhxN9cJpkHMS4yiKzSutaUIbPsjO+wx01To8ksrSPjrC6bwpONFJ063Tq+42wuiQ8HcGvr45VvjlJQ1cidM0aQWVrHR3sKzzmms0grK0l2Bl32B6SUbMkqZ05COInDAtmS5f4WpULhDM50W70EzAeqAKSU6UCfyarblqRILcx0zY+FhA3yszmKuiOuS47C1+DFx3vbO87PjO+wbYxGhQ0kOnQAadnuyXOVWVrLqu/yWDglmqevH8+U2BBWbMg+JwIsu6weIexHWllpndujH/o9MkvrKKlpYl5iBLPjw9ldcLLDicYUCk/HqZTsUsqzX0OdnsRbCJEvhDgohNgvhNjt7PVcgdVpnl5Yzay4IQhxbhdTR4QE+OpjPorbjfn4IbeKIQN97RojIQSz48PZkVvZ6VgRV2OxSJZ+cpCgAT48ec04hBA8dW0SFXXN5yR9zC6vIyY0oNOMucEBPtrcHv3Q77FZ76aakxDB7AQtkm57zvmX/FKhsOKM8SgUQswEpBDCVwjxKHoXlguYK6WcLKVMddH1nGLs0EFYXRJdGd9hi5umRnOyoaX1ISKlZEdeFdNHd2yMLokPp6HFzN7jvTu74Hs7C9h3vJqnrhnXmlp+Smwo1yVHser7PEpqTrcee7Sszu7I8rNJigricHH/y3G1ObOc8VFBDA3yJ3XEYAb6GtjiAcEQCkV3ccZ4PAA8BAwHioDJ+nafw9/HwOhw7eFoLxliZ5w95qOgqpGSmiamj+n4ejPHDMHbS/Sq36Ostonn12UxK24IN6QMb7fv8fkJWCQsX58FgNFs4Vil7ZxWtkiKDCKvssGpud7PN2oajewpOMW8RC2S3dfbi5lxYaRlVaiQXcV5izPRVpVSyp9KKYdKKSOklLdLKV0xok0CG4QQe4QQ95+9UwhxvxBitxBid0VF7z1Qp40MJXFYINGhAd06v3XMR2Y5lfXNZ/wddpzlVgL9fZgSG2p3vEdlfTNlDqR+7wp/+DyDFrOFZxdMPKdVFDM4gHtmjeKTvSc4UFRNfmUDRrPssh8oKSoIKd0bBNDbpB2twCK1LisrcxLCOVF9mtwKFbKrOD9xeCbBztKuSymXdF8cAGZJKYuFEBHARiFEppSydXZCKeUqYBVAampqr722/f668bSYnfM7LJwazcrv8li7v5j0wmrCA/0YEz6w0/MuiQ9jxYZsKuubCRngw97j1aRll5OWXcGhE7V4ewnuuWgUSy4dyyA7I7y7ysbDZXx9qJTH5icwMsy2bA/OHcOHuwv505dHuGvGSKDzSCsrs+LCGBsxiCVr9vPyrXDVxEin5D0f2JJZTmiAD5NjQlrLZuuRdFuyKojrJMRZofBEutPyeAC4CCgGdtM+HbvTKdmllMX6shz4FLjA2Wu6An8fA0H+zs19ZR3z8eHuQnbkVTGjE3+HFWvI7r3/+JGUP25k0codvJ6WR4CPN4/NT+CmqdGs+i6PS/9nC5+nF3e7K6S+2cSytYdIGBrIzy8ebfe4IH8ffnl5PLuOnWTVd7kIAWPCu9byGOTnzYcPzGBidDAPrt7LP38o6JasvUWT0cxfvs7kvZ3dk9NskWzJrmB2fHi73GXRoQHERQwirR+1wBR9i+68pkYCN6PNX24C3gc+llI67dEVQgwEvKSUdfr6FcAzzl7Xk7hpajS/W5sBnJuSxB4TooIZHTaQ8rpmrpkYyez4cGbGhRE84Iwxu2VaDMvWZrBkzT7W7DzOMz8Z32U/BGhjTn7z8QFKa5v4621T8PXu+L3i1mkxvL09n/SiGkYM6TzSqi0hAb78894LeWj1Xp767BCV9c384tKxDkex9TTZZXU8vHpv62jwqvoW/mtenENyHiiq5mRDC3MTz83cMyc+nHd2FNDYYiLA17kWo0LR2zjc8pBSVkkpX5dSzgV+BoQAGUKIO1wgz1BgqxAiHdgFfCmlXOeC63oM1jEf0Lm/w4qXl+CbX89m+xPz+MvCSVw1MbKd4QBIiQ3ls4dm8acFEzhcUstVL3/Ps18eprSmY3+IyWxhZVoul7+YRnphNX++YSJTR4R2KpO3wYvfXjMO0JIeOsoAXwMr75jKwinRvLTpKMvWZnjMXB9SSt7bWcB1r27lZEMLb909jRunDOeFjdksX5/lUMtuc2Y5XuJMN1Vb5iRE0GK28EPe+ZH8UqFoS7dfd4QQU4DFwOXA17imyyoPSHb2Op5MSIAvV08cxv7CakYM6brzvStvuwYvwe3TR3D1xEieX5fJG98f482tx7gkPpxFqTFcOi4CP+8zLYT0wmqWfnKQwyW1XJ40lKevH09UyIAuyzQnPpyH58Z1ydjYwsfgxYqbJxEW6MvKtDyqGpp58ZbJ7WR0JVJK/rnzOGlZ5UwZEcpFcWGMjwpu151Uc9rI0k8O8NXBUi4eG8YLiyYTHujH7LHh+PsY+NuWXE4bzSy7NqlLv8m3WeVMiQ09Z657gGmjtOSXW7IqmJc41KW6KhQ9jXC0f1wI8TRwLdqYjn8B66SUbom7TE1Nlbt3e8Q4Qoc43WKm2WS2+UBxJQVVDXy0p4iP9hRRUtNEaIAPC1KGc31yFJ+nF/P29nzCA/14+voJXDlhWI/K0hlvfJfHs18dITkmhJdvmWzXWd8Wq36Xjhvazhlti5pGI49/nM76jDKGBflTqkeohQT4MHPMEC6KC2dokB/L1mZQVtvEo/MTuP/i0e1yjkkpeebfh3lrWz6LL4jl2QUTbOYks1Je28QFf/6Gx+Yn8NDcOJvH3PuPHzlaXk/aY3M8rttO0XcRQuxxdhxdd4yHBcgDrKPErBcQgJRSTnJGIEc4X41Hb2O2SLbmVPLB7kI2ZpTRYrYgBNwxfQSPzk9wOhDAVXx9sITffHwAk0Wy7NokbpkWY/OBajRbeHPrMV7alE2TUYuAuz45isfmJxAz+NzW3J6CUyxZs4/yuiZ+c2Ui9140ior6ZrbnVLE1p5KtRytbjUnM4AG8cmsKKbG2W1NSSpavz+JvW3K5MWU4z980yW6G5Q9+LOTxjw/w1ZKLWyfDOpt3d+Tzu7UZbH50DqO6YDAVfQ8pJRZJuxZwT+MK49GdbqtRztxQ0fsYvLQ0J7PjwznV0MLGI2UkDA0kuZO39d7mqomRTI4N4dcfpPPEJwf5NrOcvyyc1DrCHbSutic+OciRklquSBrK41cm8tm+E7zxfR7rMkq5e9ZIHpwTR/AAHywWyarv81i+PouoEH8+fGBmawslItCfBSnDWZAyHCkluRUNHCmpZXZCeIfGVAjB41cmEuBrYMWGbBpbzKxYlGwzRHpzVjnDgvwZF2nfJzQ7PgLIYEtWOaPC3PfXamwx0dBsJjzQz20y9Ee+y67g2S+PYLJY+PCBme3quqfjcMvDk1Atj76JxSJ5c+sxlq/PIjjAhxU3J5M6IpQVG7LsdrWV1Jxm+fosPt13gpABPjw8byxp2RV8l13BNRMjeW7hRJe3sN7ceoxnvzxM7OAAXr41pZ0xbjFZmPLHjVyXHMlzN3bcGJ+3YgsxgwN4+x73RKVvy6nkF//aR2V9C8OC/JkYHUxydDCTokOYODyYUP2BZrFIGo1mGppN1DebaDKaiR8aiE8X57ZRnOFoWR3PfnWELVkVxAweQFltMxOHB/PefRfi79MzPr+2uKXbypNQxqNvc7i4lkfe30d2WT1DBvpS1dDC7dNjefzKRLuG4NCJGp798gg78qrw8/Zi2XVJ3HZBbI/5E3bmVfHL9/dTXtfczk+yPbeS297Yyao7pnLF+I79SU9/kcHqncdJ//0V5zw4ak4b2aI73W11yTmD2SJ59dujvPzNUcaED2JRajQZxbUcLKohr/LMHDJDBvrSZDTTaDRz9uNiWJA/P5s1ksXTYgkOsP2bSCnZV1jNmp3HyS6v55FLx9oMXe4PnGxo4cWN2azedZwAXwNL5o3lzpkj2HS4nIdW7+WaiZG8ujilQ1+aK1DGQxmPPk+T0czy9VnsPX6KJ68eR+rIwZ2eI6Vke24VQ4P8O00T7wpqGo088ckBvj5Uyqy4IbywaDL/930eb28vYN+yy+3O625lS1Y5P3vrR/5x97TWFCalNU38fdsxVu88Tn2zCYOX4CeTo3hwTpxLdKqsb+aRf+1na04lN6YM5083TGg31qS2ycihohoOnKihoKqBAT7eDPIzMNDPm4F+3gT6e2ORko/2FLEtp4oAXwOLUmO4Z9YoYvUowppGI5/uK2LNrkKyyuoY6Gtg8CBfCk+eZuGUaJZdm2TX4PQ1pJS8tS2fFzdpXZ0/vTCWRy6Lb9dNtTItl+e+zuQ/Zo9m6VXjelQetxoPIcQvpJQvd1bWkyjjofAUpJS8/2MhT39xGH8fL7wNXiQOC+Tdey/s9Nwmo5nkpzew+IJYbp8ey8q0PD7bfwKzRXLtpChunRbDpiPlrN5VQLPJwlUThvHgnDgmDA/ulqw786r4rzX7qDlt5JmfjGdRqu3AhK6SUVzDm1uP8UV6MWaL5IqkYQT4GfjyQAnNJguTooNZfEEs1yVH4WMQvPpNDv+blsuQgb78+YaJXJbU82HKFoskLbuCzNI6fAwCby+Bt8FLX/fC4CUwmi2YLBKT2YLRLDFZtGWT0Uxji5nGFhONLWYams2cNppIHBbE41cmdBpabrFInv4ig7d3FDA7PpynrhlncwCvlJLfrT3EP384zp8WTOD26SN66utwu/HYK6WcclbZPillijMCOYIyHgpPI6e8niVr9nG4pJbfXZvEvRd1zQl+19938UNeFc0mC/4+XtySGsN9F49u11VVVd/MW9vyeXt7PnXNJuYmhHPfxaOZMXpIl7o5apuMvLM9nxc3HSV2cACv3TbFbhRYdyitaeKdHfm8t/M4ZotkQUoUt06LtWnkDhbV8NhH6WSW1rFgchS/v258q2/FlTQ0m/h4bxFvbcvnWGX3pnMWAgJ8DAT4eRPgayDA1xtfby/SC6tJHRHK63dMJWyQ7UADs0XyxMcH+HBPET+/eBS/vXpch4baZLbw83d2k5ZdwZt3Teux7j13heouBm5Dy2/1fZtdgYBZSnmZMwI5gjIeCk+k2WRm0+FyLkuK6PKAxy/Si3nuqyPcnBrDXTNHdhh1U9tk5N0dBby59RgnG1qIDPbnJ5OHc+OU4eckqDSaLaRlVfDp/hNsOlxGs8nCNZMi+cuNEwnsoRDtZpM2J1xnureYLLy2OYfXNucQEuDDuMggWkwWjPqbv9FsocVsIdDfh6TIQJIigxgXGURiZFCnCUCLq0/z9o581uw8Tm2TickxIdx70SjmJUZgkRKTWWLUWxYmswWzReJj8MJbb4n4GLSWibeXwM/by+YD/8sDJfz6w/2EDfLjzbumkTDs3O/+l+/v598HSnjksrFdTsHT0Gxi0codHKts4IP/mNHtFmZHuMt4jEAL130OeKLNrjrgQG8OGFTGQ9GfaTKa2XC4jE/3FvHd0UrMFsn4qCBuSBlOUlQQ6w6V8kV6MacajYQG+HBdchQLUoaTEhPiUQMSM4preH5dFjWnjfh6e+Fr8MLXW3uA+xi8qKpv4UhpLdWNxtZzRgwJYEz4ILzbtLisKjW2mNmeW4WUkqsmRHLPRaO6nQWhMw4UVXPf27tpaDbx6m0prZkCmoxmHl69l01Hyvnt1Yncf8kYh65bVtvEgte2YTRLFkyOIjkmhOToEGIGD3DJb+d2h7kQYigwTd/cpWfC7TWU8VAoNCrqmvkivZhP953g4AltpkY/by8uTxrKDSnDuSQ+/LwOqZVSUlLTxJGSWg4X13KktJb8ykYsdp5fl8SHc+eMEd2ef8cRSmuauO+dH8koruXJq8ex+IJY7n93N9tyqvjjggnc0U3fRVZpHU99dpD0oprWaahDA3xaDcncxIhOMyvYw90+j5uBFcAWtNHlFwOPSSk/ckYgR1DGQ6E4l6NldeSU1zNrbJjHZA/o6zS2mPjV++msyyglPNCPqvpmlt+UzMKp0U5f22i2kFVaR3pRNemF1aQX1nC0vI6H5sbx6ysSunVNdxuPdOBya2tDCBEObJJS9lpiQ2U8FAqFp2CxSF7YmM1b246x4ubkHp3orKHZhNFs6XZ+PHelJ7HidVY3VRXOzYmuUCgU5y1eXoJH5yfwy8vjezxPVWdjh3oDZyRYJ4RYD6zRt28BvnJeJIVCoTh/6c0Eh+6k28ZDSvmYEOJGtJBdAaySUn7qMskUCoVC4bE42/bZBhjR0rLvcl4chUKhUJwPOOMwXwQsx43RVkKICqCgk8PCgMpeEMcT6c+6Q//Wvz/rDv1b/67oPkJKee7cyA5wXkdbdQUhxG5nowrOV/qz7tC/9e/PukP/1r+3dHcmOkpFWykUCkU/xdXRVl87L5JCoVAoPJ3+EG21yt0CuJH+rDv0b/37s+7Qv/XvFd1dNhmUEMIA3CqlfM8lF1QoFAqFx+Kwj0IIESSEWCqE+KsQ4gqh8TCQByxyvYgKhUKh8DS6k5J9LXAK2AFcCoQCvsAvpJT7XS6hQqFQKDwPKaVDH+Bgm3UDmiEJdPQ6vfEBrgSygBzgCXfL46DsfwfKgUNtygYDG4Gj+jJULxfAK7qeB4Apbc65Sz/+KHBXm/KpwEH9nFc48yJh8x69rHsMsBk4AmSgvZj0C/0Bf7QBt+m67k/r5aOAnbpc7wO+ermfvp2j7x/Z5lpL9fIsYH5n/wt793DD728A9gH/7oe65+v1cj+w25PrfXeU29vRtqd89AqYC4xGaxmlA0nulssB+S8BptDeeDxvrfBoE3H9t75+NVqkmwCmAzvbVIg8fRmqr1sr3i5ghn7O18BVHd2jl3WPtP4R0GaozAaS+oP+ujyD9HUf/YE2HfgAzacI8Drwn/r6g8Dr+vqtwPv6epJe5/3QHoy5+n/C7v/C3j3c8Pv/CljNGePRn3TPB8LOKvPIet8d5cxArf6pA0xt1mvd8YXbkXMGsL7N9lJgqbvlclCHkbQ3HllApL4eCWTp6yuBxWcfBywGVrYpX6mXRQKZbcpbj7N3Dzd/D2uBy/ub/kAAsBe4EG3EsLde3lq3gfXADH3dWz9OnF3frcfZ+1/o59i8Ry/rHA18A8wD/t2RXH1Nd/3e+ZxrPDyy3jvsMJdSGqSUQfonUErp3WY9yNHr9SDDgcI220V62fnMUCllCYC+jNDL7enaUXmRjfKO7uEWhBAjgRS0N/B+ob8QwiCE2I/WbbkR7W25Wp6Z4rmtvK066vtrgCE4/p0M6eAevclLwOOARd/uSK6+pjtoeQI3CCH2CCHu18s8st67Pyl8z2ErL7Jr4pI9D3u6OlruUQghBgEfA49IKWs7mLu5T+kvpTQDk4UQIcCnwDhbh+lLR3W09cLoEd+JEOJaoFxKuUcIMcdabOPQPqd7G2ZJKYuFEBHARiFEZgfHurXe9+V0IkVojlcr0UCxm2RxFWVCiEgAfWlND2NP147Ko22Ud3SPXkUI4YNmON6TUn7SiWx9Tn8AKWU1WuLR6UCIEML6stdW3lYd9f3BwEkc/04qO7hHbzELuF4IkQ/8C63r6qUO5OpLugMgpSzWl+VoLw4X4KH1vi8bjx+BsUKIUUIIXzSH2udulslZPkeLokBfrm1Tfqc+5mY6UKM3PdcDVwghQoUQocAVaH25JUCdEGK60F7n7zzrWrbu0WvoMr0JHJFSvtBmV5/XXwgRrrc4EEIMAC5DizrbDNxkQ6628t4EfCu1juvPgVuFEH5CiFHAWDRnqc3/hX6OvXv0ClLKpVLKaCnlSF2ub6WUP+1Arj6jO4AQYqAQItC6jlZfD+Gp9d4dTqHe+qBFI2Sj9Rk/6W55HJR9DVCCNl9KEXAvWt/sN2jhdN8Ag/VjBfCarudBILXNde5BC8vLAe5uU56qV8xc4K+cCdmzeY9e1v0itOb0AbSQxf36b9nn9QcmoYWpHtDlW6aXj0Z7AOYAHwJ+erm/vp2j7x/d5lpP6vploUfVdPS/sHcPN9X/OZyJtuoXuusypHMmTPvJjuqku+u9y9KTKBQKhaL/0Je7rRQKhULRQyjjoVAoFAqHUcZDoVAoFA6jjIdCoVAoHEYZD4VCoVA4jDIeCoUTCCG260uhL//Qdluh6KuoUF2FwgUIIX6FliA0AWgB0qSUG9wrlULRc6iWh0LhBEKIegCpjYQPA5YA66SUG4QQc4QQW4QQHwkhMoUQ76kWiaKv0JcTIyoUvYYQ4hG0HEmvAFcKIfzRsgOkAOPRcghtQ8vftNVdcioUrkK1PBQK1/CylPL/gAYp5ZPAJr18l5SySEppQUuzMtJdAioUrkQZD4XCBUjdeSil/EPbbaC5zWFmVGtf0UdQxkOhUCgUDqOMh0KhUCgcRoXqKhQKhcJhVMtDoVAoFA6jjIdCoVAoHEYZD4VCoVA4jDIeCoVCoXAYZTwUCoVC4TDKeCgUCoXCYZTxUCgUCoXD/D9VL+7cxEsGRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Value of pi via Monte Carlo\n",
    "n=10000\n",
    "mc_values = 50\n",
    "tic_2 = time.time()\n",
    "x_values = []\n",
    "pi_M = []\n",
    "pi_M_RMSE = []\n",
    "for i in range(1,mc_values,1):\n",
    "    M=n*i\n",
    "    if i == mc_values-1:\n",
    "        tic_2 = time.time()    \n",
    "    x_values.append([M])\n",
    "    values_call = sample_call_value(T,a,b,S0,K,r,M)\n",
    "    pi_M.append([np.mean(values_call)])\n",
    "    pi_M_RMSE.append([np.sqrt(np.var(values_call,ddof=1))/np.sqrt(n*i)])\n",
    "toc_2 = time.time()\n",
    "print( \"The last calculated value is %.2f.\" %(float(pi_M[len(pi_M)-1][0])))\n",
    "print(\"The time needed for the last run was %.5f seconds.\" %(toc_2-tic_2))  \n",
    "#The following will plot the calculated value for pi\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(x_values, pi_M, 'o')\n",
    "plt.plot(x_values, [pi_1]*len(x_values), '-')\n",
    "plt.plot(x_values, pi_1-1.96*np.array(pi_M_RMSE), '-')\n",
    "plt.plot(x_values, pi_1+1.96*np.array(pi_M_RMSE), '-')\n",
    "plt.title('i*n')\n",
    "plt.ylabel('Expected Value')\n",
    "#The following will plot the root mean square error\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(x_values, pi_M_RMSE, '-')\n",
    "plt.xlabel('i*n')\n",
    "plt.ylabel('Root Mean Square Error')\n",
    "plt.show()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2.1 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculated value with the backwards scheme was 116.5178 and it took 0.000999689102 seconds.\n",
      "The calculated value with the combinatorial scheme was 116.5178. and it took 0.000000000000 seconds.\n",
      "The calculated value with the Monte Carlo scheme was 115.6706 and it took 0.726002216339 seconds with RMSE of 2.1928.\n"
     ]
    }
   ],
   "source": [
    "print( \"The calculated value with the backwards scheme was %.4f and it took %.12f seconds.\" %(pi_1,toc_1-tic_1))\n",
    "print( \"The calculated value with the combinatorial scheme was %.4f. and it took %.12f seconds.\" %(pi_comb,toc_comb-tic_comb))\n",
    "print( \"The calculated value with the Monte Carlo scheme was %.4f and it took %.12f seconds with RMSE of %.4f.\" %(float(pi_M[len(pi_M)-1][0]),toc_2-tic_2,float(pi_M_RMSE[len(pi_M)-1][0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pricing in the N-step CRR model (5.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The computational time was 1.668977 seconds and the last approximated value is 7.012824.\n",
      "The value of the put computed with the formula is 7.012880 and it took 0.003003 seconds.\n",
      "The error between those two is 0.000056 and the time difference was 1.665974 seconds.\n"
     ]
    }
   ],
   "source": [
    "##########    Setting                                        ##########\n",
    "\n",
    "S0 = 100\n",
    "K = 90\n",
    "sigma = 0.3\n",
    "T = 1\n",
    "\n",
    "steps =np.arange(50,2001,10)\n",
    "\n",
    "#N = 10000\n",
    "#timestep = np.arange(N,N+1)\n",
    "\n",
    "##########    Price of the Put Option via CRR model          ##########\n",
    "\n",
    "#Starting the timer\n",
    "tic = time.time()\n",
    "\n",
    "#Calculating the quantities of the CRR model\n",
    "delta = T/steps\n",
    "u = np.exp(sigma*np.sqrt(delta))\n",
    "d = np.exp(-sigma*np.sqrt(delta))\n",
    "p = (1-d)/(u-d)\n",
    "q = 1-p\n",
    "n = np.size(delta)\n",
    "\n",
    "Put_value = [];\n",
    "\n",
    "for it in np.arange(0,n,1):\n",
    "    #Vector which contains all possible outcomes of the asset\n",
    "    SN = S0*(u[it]**np.arange(steps[it],-1,-1))*(d[it]**np.arange(0,steps[it]+1,1))\n",
    "    #Vector which contains all possible outputs of the put\n",
    "    P = np.maximum(K-SN,0)\n",
    "    #Using the recursion formula for pricing in the CRR model\n",
    "    for t in np.arange(steps[it],0,-1):\n",
    "        P[0:steps[it]] = p[it]*P[0:steps[it]]+q[it]*P[1:steps[it]+1]\n",
    "      \n",
    "    Put_value.append(P[0])       \n",
    "\n",
    "#Ending the timer and printing the result\n",
    "toc = time.time()\n",
    "comp_time = toc - tic\n",
    "print (\"The computational time was %f seconds and the last approximated value is %f.\"\n",
    "       %(comp_time,Put_value[n-1]))\n",
    "\n",
    "##########    Price of the Put Option via formula           ##########\n",
    "\n",
    "#Starting the timer\n",
    "tic_2 = time.time()\n",
    "\n",
    "#Incredients of the Black sholes formula\n",
    "d1 = (np.log(S0/K) + 0.5*sigma**2*T)/(sigma*np.sqrt(T));\n",
    "d2 = d1 - sigma*np.sqrt(T);\n",
    "\n",
    "#Black sholes formula for the Put\n",
    "BS_val_Put = K*norm.cdf(-d2)-S0*norm.cdf(-d1)\n",
    "\n",
    "#End timer for the formula calculation\n",
    "toc_2 = time.time()\n",
    "comp_time_2 = toc_2 - tic_2\n",
    "\n",
    "##########    Comparison outputs    ########\n",
    "\n",
    "print (\"The value of the put computed with the formula is %f and it took %f seconds.\" \n",
    "       %(BS_val_Put,comp_time_2))\n",
    "\n",
    "print (\"The error between those two is %f and the time difference was %f seconds.\" \n",
    "       %(np.abs(Put_value[n-1]-BS_val_Put),comp_time-comp_time_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD8CAYAAABkbJM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl83HWd+PHXe2YySSZJczVt0qZteoYelNKGlrOKWEoBQUEQcOVU5Ce6sv5wVdyf7nrsrrq6eK0IioIg3mgXCoiIXFKgLb1L2/ROk6Zp0+Y+Z96/P77fSSfpJJlJMpO0eT8fj3lk8pnvzLznO5N553OLqmKMMcb0xzPcARhjjDk1WMIwxhgTE0sYxhhjYmIJwxhjTEwsYRhjjImJJQxjjDExsYRhjDEmJpYwjDHGxMQShjHGmJj4hjuAoTR27FgtKSkZ7jCMMeaUsXbt2iOqWhDLsQlLGCJSCvw6omga8CVVvT/iGAG+C1wONAO3quo697YgsMk9dL+qXtXfc5aUlLBmzZohegXGGHP6E5F9sR6bsIShqtuBBW5AXuAg8GSPw1YAM93LEuBH7k+AFlVdkKj4jDHGxCdZfRiXALtUtWcmuxp4VB2rgRwRKUpSTMYYY+KQrIRxA/BElPKJwIGI3yvcMoA0EVkjIqtF5P2JDtAYY0zfEt7pLSJ+4CrgC9FujlIWXm99sqpWisg04K8isklVd0V5/DuBOwEmT548RFEbY4zpKRk1jBXAOlWtjnJbBTAp4vdioBJAVcM/dwN/A86O9uCq+qCqlqlqWUFBTB39xhhjBiAZCeNGojdHAawEbhbHuUCdqlaJSK6IpAKIyFjgAmBrEmI1xhjTi4Q2SYlIAFgGfDyi7C4AVX0AWIUzpLYcZ1jtbe5hs4Efi0gIJ6n9p6pawjDGmGGU0IShqs1Afo+yByKuK3B3lPv9HTgzkbH15rXyI+Rn+jmjcMxwPL0xxoxYtjRIhOb2Tu58dA3//fyO4Q7FGGNGHEsYEZ7eWEVTe5CahrbhDsUYY0YcSxgRfrumAoCjTe3DHIkxxow8ljBcu2saeXNvLQG/l6ONljCMMaYnSxiu1btrAXjf/Ak0tnXS2hEc5oiMMWZksYThanETxIxxmYA1SxljTE+WMFxtnU7CmJCTDsDRRuv4NsaYSJYwXK0dIQCKctIArB/DGGN6sIThausMkurzUJCZCsARq2EYY0w3ljBcbR0h0lK85Gf6AevDMMaYnixhuMI1jIDfR8Dv5YhN3jPGmG4sYbjaOkKkpjinIz/TbzUMY4zpwRKGq7UzSJrPC0B+Rqr1YRhjTA+WMFyRNYyxmX4bJWWMMT1YwnC1dgZJjahhHG2yGoYxxkSyhOFyRklF9GE0tuNs12GMMQYsYXRp6wydqGFkptIZUupbOoc5KmOMGTksYbhaO4JdNYyx7lyMI9YsZYwxXSxhuCJrGLkBJ2HU2tBaY4zpYgnD1drhTNwDyExztjpvbLMmKWOMCbOE4WrrdJYGAchKdRNGqyUMY4wJs4ThCi8NApDhJowmq2EYY0wXSxiAqtLaEbImKWOM6YMlDKA96OyFkeo2SWX4LWEYY0xPljA4sXlSuIbh9QgBv9f6MIwxJoIlDE5szxru9AanH6Op3RKGMcaEJSxhiEipiKyPuNSLyD09jhER+Z6IlIvIRhFZ2OP2MSJyUER+kKg4wVkWBE7UMMAZKdVgNQxjjOniS9QDq+p2YAGAiHiBg8CTPQ5bAcx0L0uAH7k/w74KvJSoGMPCNYzUnjUM68MwxpguyWqSugTYpar7epRfDTyqjtVAjogUAYjIImA88OdEBxfuw0iLqGFkpvqs09sYYyIkK2HcADwRpXwicCDi9wpgooh4gG8Dn01CbL3WMBrbgsl4emOMOSUkPGGIiB+4CvhttJujlCnwCWCVqh6IcnvPx79TRNaIyJqampoBxdgWpYaRleajsa1jQI9njDGno4T1YURYAaxT1eoot1UAkyJ+LwYqgfOAi0TkE0Am4BeRRlX9fM8HUNUHgQcBysrKBrSBRVtn93kYABmpXpqshmGMMV2SkTBuJHpzFMBK4JMi8iuczu46Va0CPhw+QERuBcqiJYuh0trhNkl168NIsXkYxhgTIaEJQ0QCwDLg4xFldwGo6gPAKuByoBxoBm5LZDy9CdcwIudhZKZ6aQ+G3DWmvL3d1RhjRo2EJgxVbQbye5Q9EHFdgbv7eYyfAz9PQHhdotcwwgsQWsIwxhiwmd5ARB9GRMKwFWuNMaY7SxhEXxoky1asNcaYbixhcPLig3CihmEJwxhjHJYwcGoYPo/g857ch2EJwxhjHJYwoNvmSWGZtk2rMcZ0YwkDd3vWlO4jocK77lmntzHGOCxh4CwNktajhmF9GMYY050lDKC1M3RSDcO2aTXGmO4sYQBtHcGT+jBsm1ZjjOnOEgbRaxjgdHzbNq3GGOOwhEH0GgY4CcO2aTXGGIclDJylQdKi1TDSbJtWY4wJs4SBs/hgtBpGht+2aTXGmDBLGEB758kT98AZWmubKBljjMMSBk4NI1qTVMDvpdk6vY0xBrCEATh9GNFqGE7CsBqGMcaAJQyg907vdL+XFksYxhgDWMIA4PIzC5lfnH1SeYbfR3NHEGdjQGOMGd0SukXrqeKbHzwranm630swpL3WQIwxZjSxGkYfAn4nSVizlDHGWMLoUzhhNHdYwjDGGEsYfQi4K9a22NBaY4yxhNGXcA3DJu8ZY4wljD6lh5ukrA/DGGMsYfSlq0mqw5qkjDHGEkYfAlbDMMaYLglLGCJSKiLrIy71InJPj2NERL4nIuUislFEFrrlU0RkrXu/LSJyV6Li7EtXwrA+DGOMSdzEPVXdDiwAEBEvcBB4ssdhK4CZ7mUJ8CP3ZxVwvqq2iUgmsFlEVqpqZaLijSbcJGULEBpjTPKapC4Bdqnqvh7lVwOPqmM1kCMiRararqpt7jGpSYyzG5uHYYwxJyTri/gG4Iko5ROBAxG/V7hliMgkEdno3v6NZNcuAFJ9HkRsprcxxkASEoaI+IGrgN9GuzlKmQKo6gFVnQ/MAG4RkfG9PP6dIrJGRNbU1NQMVdjhxybDb5soGWMMJKeGsQJYp6rVUW6rACZF/F4MdKtJuDWLLcBF0R5cVR9U1TJVLSsoKBiikE9I93ttWK0xxpCchHEj0ZujAFYCN7ujpc4F6lS1SkSKRSQdQERygQuA7UmI9SS2iZIxxjgSury5iASAZcDHI8ruAlDVB4BVwOVAOdAM3OYeNhv4togoTrPVf6nqpkTG2pv0FEsYxhgDCU4YqtoM5PcoeyDiugJ3R7nf88D8RMYWK9vX2xhjHDbTux8ZqT6rYRhjDJYw+pWeYvt6G2MMWMLol3V6G2OMwxJGP9L9PuvDMMYYLGH0K8NqGMYYA1jC6FfA76WlI4gzoMsYY0YvSxj9SPf7UIXWjtBwh2KMMcMqpoQhIukiUproYEairn29rR/DGDPK9ZswROR9wHrgWff3BSKyMtGBjRThhGFDa40xo10sNYx/BRYDxwFUdT1QkriQRpYTmyhZwjDGjG6xLA3Sqap1ItFWIj/9ndjXO74mqbbOIF/+0xbe2FPLmDQfv7rzPNLdxzLGmFNRLDWMzSJyE+AVkZki8n3g7wmOa8RI70oY8dUwnt5Yxa/eOsDEnHQ2VNTx/b/uTER4xhiTNLEkjE8Bc4E2nGXK64F7EhnUSBIYYMJ49PV9TCvI4NHbF3PtwmIefHk3O6obEhGiMcYkRb8JQ1WbVfWLqnqOu1HRF1W1NRnBjQQn+jBib5LaWHGc9QeOc/O5U/B4hPsuP4O0FC8Pvbw7UWEaY0zC9duHISIv4m6bGklV35OQiEaYcJNUa0fsNYwn3jxAwO/lmkXFAORnprJsznie31ZNRzBEitemvxhjTj2xdHrfG3E9DbgWGDWTEgIp8TdJrdt3jHOn5TMmLaWrbMW8Qp58+yCv7zrK0llDv5WsMcYkWr8JQ1XX9ih6TUReSlA8I064htESYw2jtSNIeU0jl84d36186awCMvxentlcZQnDGHNKimXiXl7EZayILAcKkxDbiJDq8yAS+8S9dw41EAwpcydkdytPS/FyyezxPLelms6gLTNijDn1xNIktRanD0NwmqL2AHckMqiRRETi2kRp88E6AOZOGHPSbZfMHsfKDZVsrapnfnHOkMZpjDGJFkuT1NRkBDKSBfxemmNsktpSWU92egrFuekn3bZkqrO9+Zt7ai1hGGNOOb0mDBG5pq87quofhj6ckSktxUtrjDWMLZV1zJ0whmgz4wuz05icF+DNPbV89KJpQx2mMcYkVF81jPf1cZsCoyZhxLpNa0cwxDuHGrj1/JJejzmnJI8Xtx9GVaMmFWOMGal6TRiqelsyAxnJ0lO8MY2S2lXTSHtnKGr/Rdjiqbn8fl0Fu2oamTEuayjDNMaYhIql0xsRuQJneZC0cJmqfiVRQY006f7YOr33HW0GYNrYzF6PWdzVj3HMEoYx5pQSy7DaB4AP4awpJcB1wJQExzWixFrDqDreAkBRTlqvx5TkBxibmcqavbVDFp8xxiRDLGtUnK+qNwPHVPXfgPOASYkNa2QJ+H0xrSVVWdeK3+chP8Pf6zEiwlnF2Wxyh98aY8ypIpaE0eL+bBaRCUAH0O9QWxEpFZH1EZd6EbmnxzEiIt8TkXIR2SgiC93yBSLyuohsccs/FO8LG0ppKd6Y9vSuPN7ChOy0fjuz507MZldNY9x7bBhjzHCKpQ/jKRHJAb4FrMMZIfVQf3dS1e3AAgAR8QIHgSd7HLYCmOlelgA/cn82Azer6k43Sa0VkedU9XhMr2qIOaOkYqhhHG+hKPvk+Rc9nTkxm5DCtqoGFk3JHYoQjTEm4WJZ3vyrqnpcVX+P03dxhqp+Kc7nuQTYpar7epRfDTyqjtVAjogUqeoOVd3pPn8lcBgYtgWY0v0x9mHUtTIhp/+EMW+iM4pqS6U1SxljTh2xdHpvEJH7RGS6qrap6kC+5W7A2Xypp4nAgYjfK9yyyOdfDPiBXQN43iGR7jZJhUInrfLepTMYorq+lQl9dHiHFY5JIz/Dz6YKSxjGmFNHLH0YV+GsIfUbEXlLRO4VkcmxPoGI+N3H+G20m6OUdX0ri0gR8AvgNlWN2okgIneKyBoRWVNTUxNrWHHp2hOjs/daRnVDGyElpiYpEWHuxGw2V9YPWYzGGJNosTRJ7VPVb6rqIuAmYD7OAoSxWgGsU9XqKLdV0H3EVTFQCSAiY4CngX9xm6t6i+9BdyfAsoKCxLRapcewJ0Z4SG0sNQyAeRPGsLO6Ia6NmWJ1qK6VP759kOe2HOJIY9uQP74xZnSKdeJeCXA9znyMIPDPcTzHjURvjgJYCXxSRH6F09ldp6pVbq3kSZz+jWg1k6Tq2hOjj4RxsCth9F/DAJg3MZvOkLKjumHIFiJs7wzxlae28Njq/V1lXo/wsYumcc97Z5LmJj5jjBmIWLZofQNIAX4DXKeqMW9MLSIBYBnw8YiyuwBU9QFgFXA5UI4zMiq8HMn1wFIgX0RudctuVdX1sT73UArXMPqqDVTVOducF2XHVsMoLXRmeW8/NDQJo60zyC0Pv8nq3bXcen4JH1xUTFtniF+/tZ8HXtrF67uP8tgdi8mK2AXQGGPiEUsN4xZVfWcgD66qzUB+j7IHIq4rcHeU+z0GPDaQ50yEgL//JqnK4y1kpfli/kKekhfA7/Ow83DjkMT4nT/vYPXuWv7rurP4oLuXOMCiKbm854xxfPKXb3PHz9fwyO2Lu2pMxhgTj1j6MAaULE4n4RpGX0NrK4+3MiGGDu8wn9fDjIJMth9qGHR8r+86yoOv7OamJZO7JYuwy+YVcf8NC3hrXy1feWrroJ/PGDM6xTJKatSLpQ+jur6Vwhibo8JKC7PYUT24hKGq/PuqbRTnpvMvV8zu9bgr50/g40un88Sb+3l286FBPacxZnSKZR5Gaixlp7OuhNFHDaO2qb3PNaSimTU+i6q6VupaOgYc26vlR9h0sI673z2DgL/vFsbPLJvF/OJs7ntyE3XNA39OY8zoFEsN4/UYy05bgRTni7ivPoxjze3kBOJLGKWFzjLoOwdRy/jhi+WMH5PKBxZO7PdYv8/Df1xzJseb27n/hR0Dfk5jzOjUa8IQkUIRWQSki8jZIrLQvbwbCCQtwhEgze+cpt5qGK0dQZrbg+RlxDcCadZ4d6TUABPGtqp6Vu+u5WMXTSPVF1tH9twJ2dy4eDKPvr5v0M1hxpjRpa82jOXArTiT6b4TUd4A3JfAmEaccFNPSy8LEB53m3firWFMzEknw+9lxwA7vv+0vhKfR7hm4ckd3X2599JSVm6o5Nt/3s6PP1I2oOceqOr6Vl4rP8L+2mYaWjvJDaQwrSCTJVPzyM8cVS2dxpxy+tqi9RHgERG51l14cNTqGiXVHn2J82PN7QDkxdmHISLMHJ/Fjur4h9aqKv+7oZILZ46N+3lzM/zcfsFUvvvCTrZU1jF3Qnbczx+vt/bW8oO/lvPSjhPLt0TulS4CF80s4Nbzp3Bx6Tjb79yYESiWeRjzRGRuz8LRtEWr1yP4fR6aO6LXMMIJIycQ/6S4meMyu32Jxmrd/mMcPN7C/710Vtz3Bbj9wqk8/Noe/vv5nfzklsTVMpraOvn6qm388o39jM30c897Z7Jsznhmjc8ixeuhtSPI1qp6/vbOYX67toLbf76GxSV5/OtVc5nTx97oxpjki6XTuxFoci9BnLWhShIY04iUnuKltZdO72NNTpNUvP/pA0wryORwQxv1rfGNWlq5vpJUn4dL5xbG/ZwA2ekp3HHhVP6yrTphfRk1DW1c/+PXeeLN/dy5dBqvfu493PPeWcydkE2K1/nopaV4WTg5l89cWsor/3wxX//APHYfaeTqH77KD18s73OFYGNMcsUyce/bEZevA++mxxLko0Fk80lP4RpGbpx9GADTCzIA2F3TFNf9Xtxew0UzC8hMjWk5sKhuPq+EVJ+Hh1+NZy3J2Byub+WDD/yd3TVNPHzrOdx3+ex+17LyeT18eMkUnv+nd3Hp3EK+9dx27vzFWhriTKbJoqq0tAc53NBKTUMbTW2dOIsXGHN6Gsi3TQCYNtSBjHTpKb1vonSsaeBNUtPHOUNrdx1uZMGk2NaU2n+0mf21zdxxYb875fYpL8PPtYuK+d3aCu5dXsrYIep0bmrr5PZH3qKmoY3HP7aEhZPj21UwN8PPD248m3Om5PLVp7dx/Y9X8+jtiynIGt5O8ca2Tl585zCrdx9l08E6dtc00djWvZky4Pcyc1wmcyZkM3fCGJbOLGBy/qgaVGhOY7EsPriJE3tUeHF2vhs1/Rdh6X5vrzO9jzV3kOH3xjy0NdLkvAA+j7CrJvaO71fLjwBwwYyxcT9fT3dcOJVfvrGfx1bv4573Dqw/JJKq8tnfbWBrZT0/uaUs7mQRJiLcesFUphZkctcv1nLdA3/nV3eeF/ds+qGwdl8tP//7Pp7bcoj2zhCZqT7OmpTNBxcVM35MGplpPlClqT1IdX0r71Q18PTGSp5401k1eE7RGC6bV8g1CydSnGvJw5y6YqlhXBlxvROoVtX+N7g+zfRVwzje3E7uAPovAFK8HqbkB+JMGDUUZad1NWcNxvSCTJbOKuBXbx7gkxfPwOcd3Goxv11bwapNh/jcZWfwnjPGDzq+d80q4LGPLuHmn77BzQ+/wW8+fl7cw5cHalNFHd949h1eLT9CVpqPmxZP5or5RZw9Kaff86Sq7D3azF+2VvPslkN85/kd3P+XHVw2r5A7Lpx2yuzl3toRZFdNI9X1TrNbTUMbbZ0hgiElpJCW4iE7PYUxaSmMH5PGlPwARdlpg/4cmZGp14QhImnAXcAMYBPw09GYKMLS/V4aWqO//Nrm9gH1X4RNK8iMuQ8jGFJeKz/Ksjnjh2zo6T8smcydv1jLX985POBOdIADtc3828otnDstjzuXDl2r5aIpuTx0Sxm3PvwWt/38LR7/6JJ+l0EZjOb2Tr757HYefX0vuQE//3LFbG5aMjmu5xQRpo7N4GNLp/GxpdM4eLyFR1/fyxNv7GfVpkMsmZrHfZfP5qwYmyGT5eDxFl7eUcPru46ytaqePUeaCEYZeODzCB4R2oMnDzX3eZzXfmZxNvMnZnP25FzmTczG67Gh0qe6vv4CHgE6gFdwRkbNAT6djKBGovQULzUN0XevO9bcMaD+i7DpBZn8bfthOoOhfv8z21pZT11LBxcOQXNU2HvOGEfhmDQef2P/oBLG157eSkjh29cvGPIvh/Onj+X7N53N/3lsLR//xVp+ess5+H1D/1/s1sp6PvXEOnYfaeIj507h3uWljBmCPUQm5qTzhRWz+cf3zOQ3aw7wg7+Wc/UPX+Oqsybwz5eVDmtT1f6jzfxuXQVPb6xkl/uPy/gxqZw5MYfL5xVSWjiGCTlpFGSlMjYztdvghc5giIbWTupaOqiqa2V/bRN7jzaz41ADL+84wh/WHQRgTJqP86ePZemsApbNGT/s/VHxau0Isrumif21TeyvbabyeCtHGts42tjOseZ22jtDtLmX9s4gIkJ6ipd0v5e0FC/pKR7yMvyMG5PG+Kw0xo1JpXBMGtMKMijODZwyybSvhDFHVc8EEJGfAm8mJ6SRqc9RUk3tlAyiY3N6QQYdQeXAsRamju27mWnd/mMALJ6aN+Dn68nn9XDD4kl894WdVBxrHtCX1ys7a3huSzWfXV7KxBh3HYzX8rmF/Oe18/nn323kq09t5avvnzekj//s5iru+fV6xqSl8PgdSzh/CJNyWEaqj9sumMoHFxXz45d285NXd/P8Vue83XJ+SdK+OIIh5ZnNVTz6+j7e3FOLCJw/PZ8bF0/mXbMKmDEuM6YarM/rITfDT26Gn5KxGZw3/cT2N6rKofpW3tp7jNd2HuHV8iM8u+UQX/zjJhaX5HHF/CLeN3/CgJtzEyUYUrZV1fPmnlo2VBxna2U9u2oaiaxoZaX5KMhMJT/Tz6S8AKk+D36fx/np9aA4Saa1I0RLR5CW9iAHj7fy9v7jHHUHyYT5fR6m5mcwfVwGs8ZncdakHM4qzhnQMP1E6ythdI1lVNXO0T7zNt3fxyipQTZJRY6U6i9hbDhwnIKs1Jh39ovVtQuLuf8vO/nT+kruvnhGXPcNhZSvPbWNKfmBQY/c6s/1ZZPYWd3AQ6/sYcGkHK6Nsv/HQDzy9718eeUWFkzK4aGbyxL+H3BWWgr3Li/lxiWT+ZcnN/GVp7by1MZKvvnB+cwYl5Ww5+0Ihli5vpIf/q2c3TVNTMkP8NnlpVyzcCJFceznEgsRoSg7navOSueqsyagquyobmTVpipWbariS3/awtee2sayOeO5rqyYpTML8AzTf9qH61t54Z3DvLDtMG/sPkqDO/qtcEwacyc4gxZmjc+iJD+DyfkBstMHXuts7wxxpLGNyuMt7K5pYldNI7tqGtlaWc8zmw8RHpk9OS/A/OJsFk/N4/zpY5lekDHsKyD0lTDOEpF697rgLEJY715XVR1V03DTU3xRR0l1uFXywSSMqflOkth7tP9+jPUHjrNgUs6Qf3Am5QVYXJLH79dV8Il3T4/r8Z/ZfIjt1Q1894YFSdk3/HOXncGmg3Xc9+QmSguzmDdxcEubPPTybr6+yvni+v6NZyd17/OJOek8fOs5/HH9Qf7tf7dyxfde5StXz+X6sklD/h6/Vn6EL6/cQvnhRs4ozOIHN53NinlFSavViAilhVmUFmbxT8tmsa2qnt+sOcAf3z7I05uqmDY2g9svdGpfyXgPjjS28b8bKvnj+ko2HDgOOO/HlWcVsWRqPoun5jEhAbVlv8/DhJx0JuSkU1bSvaWgsa2TTRV1bKg4zsaK46zdd4ynNlYBTjPh+dPHcv70fN5dOm5YmvX6WkvK9vGMkO730NIRRFW7/SGHFx7MjXOl2kg5gRTGpPn6TRh1zR3sPtI0ZP9V9/SBhRP5wh82sbGiLubO2FBI+e4LO5hekMGV8yckJK6efF4PP7hpIVd+71XuemwtT33qwgGPnHps9T6+vmobV8wv4v4PLeiagZ5MIsIHzi7mghlj+cyvN/C532/i9V1H+doHzhzUxMywqroWvvbUNp7eVMXkvAA//sgiLh3CQRMDNbtoDF9+31w+v+IMnt18iJ++uod/+eNmvvP8Dv7h3CncfN6UIZsbFBYKKX/bcZjHVu/npR01BEPK3Alj+OzyUi6ZPY7S8VnDel4yU32cNz2/q2lPVdlf28zfdx3ltfIjvLyjhiffPogInD0ph2VzClk2ZxzTC2JrQhwsOZ1mppaVlemaNWvivt8998D69X0fc/B4Cwdqm1kyNa/bG9PS3smGijpmjssc1Gqrmw/W4fUIs4t6r7jVtbSzraqB2UVjBlUl7k1nKMS6fccYNyaNkvzYhuzWNrWxo7qRGeMyh/yPuz+NbR1sqawnL+Bn5vhMnMpv7MKx5wZSmDXMXxRhilJ5rIWKYy2kpXgpLcwa1H/bRxvb2H2kCVVlYk46RTnpeEbA64xOqW/tpOp4C8eaO/AIFGanMSE7fdDDdIMh5XBDK9X1rbR2hPB7PYzNSqUg0096AkfcDT2lqS3IsWans72pzWn1GFfSwu4XJwzoPInIWlWNaUG5U+lMDatwrT2oii/iD67D7Qkb7Ac6LcXT1W7am0Z3WG9mamIqfz6Ph9yAn6ONbUzJC8T0BXqorpXUFA/5mcnvoMtMTWFSboD9tc3kNrbHlbCa2jspP9xEZpqPmSMkWQAIwsTcAFlpKew83MDmg3XMHJ9Jdnp85zcYUvYebaKmoc15jQWZpCaxqW1ghDFpKYwpTKGlI8jBYy1UHm+lur6Nouw0irLT8Hri+zsLhkIcqm+lqq6VzqCSleZjUm6AvAz/iHnP4yNkpPrISPVRnBugvTPIseYO8vMkKXNfLGEA99/f/zG/fKOG+57cxO++cEm32cbPbj7KXY+t5bF/vHBQy4R/58+V/ODFcv781RW9Dhf96CNb2HOkiRf+77sH/Dz9+cvWZj766Bq+dEsZl8zue+Ld5oN1XPn9V/n3K2bz0YuGZyJaMJTGh378NtsgxIsqAAAVD0lEQVSrG3j8nqUxjdCqa+ng8u++wlmq/OmTFzAuayR+caRwoDaNjz6yhvKaRu65cg63nF8S0z33HGnijp+/RcrRJv794hn84yUzT8GJdF4gk+2HlO88v53ntlTTGEjh7otncPN5Jf0OqW7tCPLzv+/lf14sp6O1k+vPGMen3jODsycnfin/5PK6l+SsgHCqfYqGTcDd17u5xyZK4VVmBztWf0p+BiGFimPNvR6zpbKeMwfZwdufd5UWkJfh5w9vH+z32Ef+vpf0FC/XlU1KaEx98XqE71y/gFBIufc3G2Ja3fZrT23lUH0rP/qHRYzLSv5SI7GalBfg9584n4tLC/jyyi38+6pt/b6+1buP8oH/eY3jLR388mPn8plLS0/BZHFCaWEWP/5IGX+6+wLmTczma09vY/n9L/PnLYeiLvQY3ifmvd95if985h3KSvJ46lMX8vCt53D2AJepMSecup+kJAu3I/ccWhtuJspKG1xlrWSsM/dh39HoCSM8Maq0MLGD01K8Ht43v4jnt1b3ueR6XUsHKzdUcs3CiQnpT4nH5PwAX37fXF7ffZRHXt/b57F/2+7su/HxpdNiXuxxOGWm+vjxR8q4+bwpPPjybu793QY6osyuBvj92go+8tM3yM/w8+QnzufcaflRjzsVnTUph1/csYSf3XYOXo9w5y/WctNDb7C1sr7rmD1HmrjpoTf41BNvk5nq4/GPLuHhW88Z9Cg6c4I1ScUoXMPoObS2ye13yBjkaJYpbifzniNNXBzl9p3unhWlhZmDep5YXH32RB55fR9/2Vrd6/avT2+soq0zxIfOGb7aRaTryopZtbmK/3puO8vnFkYdDlnf2sEX/rCJmeMy+fR7Zw5DlAPj9Qj/dtVcCjJT+fbzOzjW1M7/fHgR6f4TfRLheSQXzMjnfz68aNiTeKJcXDqOC2eM5Zdv7Oe//7KDK77/CtcvmkRxbjrff7GcVJ+Hr71/HjcunnzKzJ4+lVgNI0bhP86TahhtnaT6PIMejpmf4Scz1ce+XobWbncTxqzxiZvUFbagOIei7DRWbTrU6zG/W3uAWeMzE95EFisR4atXzyOk8KU/bY7aXPEfq7ZRXd/Kt647a0ArCw8nEeFTl8zkP645k5d21HDLz97sWlr94Vf38OWVW1g+dzw/u3XxaZsswlK8Hm45v4SX7r2Y2y+Yyh/eruDbz+/gvbPH8cJn3sU/nDvFkkWCJKyGISKlwK8jiqYBX1LV+yOOEeC7wOVAM3Crqq5zb3sWOBd4VVUjV8wdFuF9vXsuD9LY1jkkY+VFhJKxAfb20iS141ADGX5vwpbdiOTxCMvnFvLLN/dHfX27ahpZt/84911+xogaaTIpL8Bnls3i66u28ezmQ6w4s6jrtld3HuGJNw+cMk1Rvblx8WQyUn3806/Xc9X3X+W86fk8/sZ+Vswr5Hs3nj0s80iGS3Yghf935RxuOa+E2ub2U/p9PVUk7NOlqttVdYGqLgAW4SSEJ3sctgKY6V7uBH4Ucdu3gI8kKr54hWsYrVFqGJmD7L8Im5Kf0WcNY1Zh8oZ/Xn5mEe2dIV585/BJt/3x7YN4BN6/YORtvHjbBSXMnTCGL6/c0tUH09jWyed+v5FpYzP4p2WD3/NjuF111gR+dquz+OLjb+znijOLRl2yiDQ5P2DJIkmS9Qm7BNilqvt6lF8NPKqO1UCOiBQBqOoLQGI2mx6AE6OkTu7DyBiiiT8l+QEqjrWc1Kmpqmw/1EBpEpqjwhZNyWVsZirPbu7eLKWqPL2pinOn5TNuzMgbYeTzeviPa87kSGMb//XcdgC+8cw7VNa18K3r5id12Y9EWjqrgGc+fRHPfPqiUZ0sTHIl61N2A/BElPKJwIGI3ysYofuFh5ukenZ6N7QObQ2jM6QcPNbSrbymsY1jzR1J6b8I83qE5XPH89d3Dnd7zTsPN7K7pqlbc89IM784h4+cO4XHVu/j4Vf38IvV+7jt/KksmjJ0K/yOBCLOygDWXm+SJeEJQ0T8wFXAb6PdHKUsrrVKROROEVkjImtqamoGEmJMeuv0bmofmj4MoGs5jp5rSu045OzGV1qYvIQBsGJeES0dQV7aceK8rtpUhQgsnzv43fQS6TPLSskN+PnKU1u7VmQ1xgxOMmoYK4B1qlod5bYKIHJcZjFQGc+Dq+qDqlqmqmUFBQWDCLNvfq8Hj5xcw2hsHcKE0ctcjPLDTsvczHGJH1Ibacm0PHIDKTy7uaqr7JlNhzinJG9ET3gDp0P0i1fMJi3Fwzevnd9tCKoxZmCSkTBuJHpzFMBK4GZxnAvUqWpVL8cOKxEh4PdFGSUVHPQcjLCCzFQCfu9JNYy9R5vJ8HuTvpxxitfDsjnjeWHbYdo6g+w72sT26gYuG8SufMl0zcJi1n/pUpacRhPYjBlOCU0YIhIAlgF/iCi7S0Tucn9dBewGyoGHgE9EHPcKTjPWJSJSISLLExlrLNJSTt5EqbGtY9CzvMNExB0p1b2Gse9oE1Pyh2fzlBXzimho6+SVHUe6Rky954xxSY9joE6XTm5jRoKEzvRW1WYgv0fZAxHXFbi7l/telMjYBiLg99ISsZZUZzBEa0doyEZJgTNSKjxJL2zf0WbOKEpu/0XYhTPHMjYzlV+9dYCOYIhpYzMo6WdXQGPM6cnG4sUhvUcNI7wW/VCNkgJnpNSB2mY63aG1ncEQ+2ubu5YOSbYUr4fryor56zvVvL7rKBefQrULY8zQsoQRh3S/t1sfRmP70O9PUZIfoCOoVNW1AlB5vJXOkFKSHxiy54jXjedMJqTQHgxxcaklDGNGK0sYcUhP8Xab6X1iQ6OhW7sn3NwT7vgO/4x1B7xEmJwfYOmsAjL8Xs6ZaktEGzNa2Wq1cQj4vRyqP7Hkd2PXSrVDWcNwE8aRJi6aWdC1VMhw9xt849ozqa5vO+UW7TPGDB1LGHFI83fvwwgnjKEaJQUwLiuVVJ+HA+5s7z1HmklP8TIuyUNqeyrKTqcoO/ELHxpjRi5rkopDIMVLa3tkp/fQ7IURyeMRinPT2e8OrXWG1Ma2v7YxxiSSJYw4pPu9NEftwxjaitrkvAD7a52Esfdo07D2XxhjTJgljDik+73dlgYJN0klImEcqG2mvdMZUju1wBKGMWb4WcKIQ3qKl7bOEMGQsz5iYwKapMDZCKihrZM1e2vpCCpzihK7j7cxxsTCEkYcAj02UWoaou1Ze5qc58y5eHaLsxfFnAmWMIwxw88SRhx6btPaMETbs/Y02Z2k99yWQ6SneK0PwxgzIljCiEO6u2ZUZA1jKJcFCZuU6ySM6vo2zijKsg1yjDEjgiWMOPSsYTS2Dt32rJEyUn2MzfQDWP+FMWbEsIQRh3S/c7rCk/caE1TDAKfjG2C2JQxjzAhhCSMO6SlOcmh2Fx1sTFAfBpzo+LYOb2PMSGEJIw7pUUZJJSphzCjIJNXn4Ywk7+NtjDG9sbWk4hAeVtvVh9HWOeRzMMLuuGgql80rJJCAPhJjjBkIq2HEIdzpHZ7t3dIe7EoiQy3g9zFzvNUujDEjhyWMOISbpFo6gqgqrZ0h0lLsFBpjRgf7totDZA2jI6gEQ9pVZowxpztLGHGInIcRHlqbZgnDGDNKWMKIg8cjpPo8tHYEabOEYYwZZSxhxCng93arYViTlDFmtLCEEaf0FGeb1q6EkaBRUsYYM9JYwohTuruvd2tHCMBGSRljRg37totTeNe98FwM68MwxowWCUsYIlIqIusjLvUick+PY0REvici5SKyUUQWRtx2i4jsdC+3JCrOeAVSfLS0B7uWB7E+DGPMaJGwdSdUdTuwAEBEvMBB4Mkeh60AZrqXJcCPgCUikgd8GSgDFFgrIitV9Vii4o1Vmt9LXUtHV8KwGoYxZrRIVpPUJcAuVd3Xo/xq4FF1rAZyRKQIWA48r6q1bpJ4HrgsSbH2KZDipdVGSRljRqFkJYwbgCeilE8EDkT8XuGW9VY+7NL9Xpo7Ors6vW2UlDFmtEh4whARP3AV8NtoN0cp0z7Koz3+nSKyRkTW1NTUDDzQGDmd3qETM719ljCMMaNDMmoYK4B1qlod5bYKYFLE78VAZR/lJ1HVB1W1TFXLCgoKhijk3qWneGlp7zzRh+G3gWbGmNEhGd92NxK9OQpgJXCzO1rqXKBOVauA54BLRSRXRHKBS92yYRdw52G0tAfxCPi9ljCMMaNDQnfnEZEAsAz4eETZXQCq+gCwCrgcKAeagdvc22pF5KvAW+7dvqKqtYmMNVZpKV5CCnUtHaSleBGJ1npmjDGnn4QmDFVtBvJ7lD0QcV2Bu3u578PAw4mMbyDCGybVNrfbCCljzKhi7SlxCieJY03tNgfDGDOqWMKIU3gYbW1Tu60jZYwZVewbL05dNYzmdpuDYYwZVSxhxCmcJI41ddgcDGPMqGIJI07hTu/2YMhqGMaYUcUSRpwiO7qt09sYM5pYwohTwH9iJLIlDGPMaGIJI06Rcy/SbZSUMWYUsW+8OEX2W9jEPWPMaGIJI07p1odhjBmlLGHEye/z4PM460dZwjDGjCaWMAYgXMuwYbXGmNHEEsYAhBNFms9OnzFm9LBvvAEIJwyrYRhjRhNLGAMQbpKyPgxjzGhiCWMAupqkLGEYY0YRSxgDEF5PyuZhGGNGE0sYA2CjpIwxo5EljAFId9eTsuXNjTGjiSWMAQivIZXut9NnjBk97BtvAMIr1lqntzFmNLGEMQBpNqzWGDMKWcIYgIANqzXGjEK+/g8xPV0xvwivR8hMtdNnjBk9rIYxANMLMrn74hnDHYYxxiSVJQxjjDExSWjCEJEcEfmdiLwjIttE5Lwet+eKyJMislFE3hSReRG3fVpENovIFhG5J5FxGmOM6V+iaxjfBZ5V1TOAs4BtPW6/D1ivqvOBm93jcRPHx4DF7v2uFJGZCY7VGGNMHxKWMERkDLAU+CmAqrar6vEeh80BXnBvfwcoEZHxwGxgtao2q2on8BLwgUTFaowxpn+JrGFMA2qAn4nI2yLyExHJ6HHMBuAaABFZDEwBioHNwFIRyReRAHA5MCmBsRpjjOlHIhOGD1gI/EhVzwaagM/3OOY/gVwRWQ98Cngb6FTVbcA3gOeBZ3ESS2e0JxGRO0VkjYisqampScwrMcYYk9CEUQFUqOob7u+/w0kgXVS1XlVvU9UFOH0YBcAe97afqupCVV0K1AI7oz2Jqj6oqmWqWlZQUJCo12KMMaNewhKGqh4CDohIqVt0CbA18hh3FJXf/fWjwMuqWu/eNs79ORmn2eqJRMVqjDGmf6KqiXtwkQXATwA/sBu4DfgQgKo+4A6zfRQI4iSTO1T1mHvfV4B8oAP4jKq+EMPz1QD7erl5LHBkUC8ocSy2gbHYBsZiG5jTNbYpqhpT80xCE8ZIIiJrVLVsuOOIxmIbGIttYCy2gbHYbKa3McaYGFnCMMYYE5PRlDAeHO4A+mCxDYzFNjAW28CM+thGTR+GMcaYwRlNNQxjjDGDcNonDBG5TES2i0i5iPScaZ6M558kIi+6q/VuEZFPu+X/KiIHRWS9e7k84j5fcOPdLiLLExzfXhHZ5Mawxi3LE5HnRWSn+zPXLRcR+Z4b20YRWdj3ow8qrtKIc7NeROpF5J7hPG8i8rCIHBaRzRFlcZ8rEbnFPX6niNySwNi+5a4UvdFdFTrHLS8RkZaIc/hAxH0WuZ+Hcjd+SVBscb+Pifhb7iW2X0fEtVeclSiSet76+N4Y3s+bqp62F8AL7MJZ18qPs8TInCTHUAQsdK9nATtwFl38V+DeKMfPceNMBaa68XsTGN9eYGyPsm8Cn3evfx74hnv9cuAZQIBzgTeS+D4ewllrbNjOG85imguBzQM9V0AezpykPCDXvZ6boNguBXzu9W9ExFYSeVyPx3kTOM+N+xlgRYJii+t9TNTfcrTYetz+beBLyT5vfXxvDOvn7XSvYSwGylV1t6q2A78Crk5mAKpaparr3OsNOEu8T+zjLlcDv1LVNlXdA5TjvI5kuhp4xL3+CPD+iPJH1bEayBGRoiTEcwmwS1V7m5QZji2h501VX8ZZpqbn88ZzrpYDz6tqrTqTVJ8HLktEbKr6Z3VWewZYjbOwZ6/c+Mao6uvqfNs8GvF6hjS2PvT2Pibkb7mv2NxawvX0s8pEIs5bH98bw/p5O90TxkTgQMTvFfT9ZZ1QIlICnA2E19f6pFt9fDhctST5MSvwZxFZKyJ3umXjVbUKnA8uMG6YYgu7ge5/tCPhvIXFe66GK87bcf4DDZsqzirSL4nIRW7ZRDeeZMUWz/s4HOftIqBaVSPXsUv6eevxvTGsn7fTPWFEa0cclmFhIpIJ/B64R531sn4ETAcWAFU4VV9IfswXqOpCYAVwt4gs7ePYpJ9PcdYauwr4rVs0Us5bf3qLZzjO4RdxVnt+3C2qAiars4r0Z4BfirN/TTJji/d9HI7390a6/6OS9PMW5Xuj10N7iWFIYzvdE0YF3ffRKAYqkx2EiKTgvOmPq+ofAFS1WlWDqhoCHuJE80lSY1bVSvfnYeBJN47qcFOT+/PwcMTmWgGsU9VqN84Rcd4ixHuukhqn28l5JfBht7kEt7nnqHt9LU7fwCw3tshmq4TFNoD3MdnnzYez6OmvI2JO6nmL9r3BMH/eTveE8RYwU0Smuv+p3gCsTGYAbjvoT4FtqvqdiPLItv8P4GwaBU58N4hIqohMBWbidKglIrYMEckKX8fpJN3sxhAeTXEL8KeI2G52R2ScC9SFq8cJ1O2/vJFw3nqI91w9B1wqzn72uTjn/LlEBCYilwGfA65S1eaI8gIR8brXp+Gcq91ufA0icq77ub054vUMdWzxvo/J/lt+L/COqnY1NSXzvPX2vcFwf94G05N/KlxwRg/swPlv4IvD8PwX4lQBNwLr3cvlwC+ATW75SqAo4j5fdOPdzhCMUukjtmk4o002AFvC5wdnleAXcPYgeQHIc8sF+KEb2yagLMHnLgAcBbIjyobtvOEkriqcFZQrgDsGcq5w+hPK3cttCYytHKf9Ovy5e8A99lr3/d4ArAPeF/E4ZThf3ruAH+BO7k1AbHG/j4n4W44Wm1v+c+CuHscm7bzR+/fGsH7ebKa3McaYmJzuTVLGGGOGiCUMY4wxMbGEYYwxJiaWMIwxxsTEEoYxxpiYWMIwxhgTE0sYxhhjYmIJwxhjTEz+P/fRsb7pzmIHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.plot(steps, Put_value, steps, [BS_val_Put]*len(Put_value), 'b-')\n",
    "plt.ylabel('Put value')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5.2.2 (4 points)\n",
    "\n",
    "See explanation above.\n",
    "\n",
    "Note that you will have to specify a computationaly reasonable backwards scheme and a suitable Monte Carlo sampling routine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2.2 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The payoff of the knock-out option is given by\n",
    "$$\n",
    "    \\xi^{\\text{k.o.}}=\\frac{(S^1(T)-K)_+}{(1+r)^T}\\mathbb{1}_{\\{\\sup_{v\\in \\{0,\\dots,T\\}}S^1(v)\\leq B\\}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2.2 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "S0 = 120\n",
    "a = -0.5\n",
    "b = 0.5\n",
    "r = 0.1\n",
    "B = 200\n",
    "K = 120\n",
    "T = 300\n",
    "#Unique equivalent martingale measure:\n",
    "p_star = (r-a)/(b-a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backward formula for arbitrary $T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculated value is 0.0000.\n",
      "The time needed was 0.0160 seconds.\n"
     ]
    }
   ],
   "source": [
    "#backward formula of the expected value for general T\n",
    "tic_1 = time.time()\n",
    "#start the time\n",
    "asset_val = S0*((1+b)**np.arange(T,-1,-1))*((1+a)**np.arange(0,T+1,1))\n",
    "#Vector which contains all possible outcomes of the asset\n",
    "asset_val[asset_val>B] = 0\n",
    "call_val = np.maximum(asset_val-K,0)\n",
    "#Vector which contains all possible outputs of the knock-out call\n",
    "#Next we use the recursion formula for pricing in the CRR model:\n",
    "for t in np.arange(T-1,-1,-1):\n",
    "    asset_val_temp = S0*((1+b)**np.arange(t,-1,-1))*((1+a)**np.arange(0,t+1,1))\n",
    "    asset_val_temp[asset_val_temp<=B] = 1\n",
    "    asset_val_temp[asset_val_temp>B] = 0    \n",
    "    call_val[0:t+1] = np.multiply(p_star*call_val[0:t+1]/(1+r)+(1-p_star)/(1+r)*call_val[1:t+2],asset_val_temp)\n",
    "pi_1 = call_val[0]\n",
    "toc_1 = time.time()\n",
    "#Ending the timer\n",
    "print( \"The calculated value is %.4f.\" %(pi_1))\n",
    "print(\"The time needed was %.4f seconds.\" %(toc_1-tic_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combinatorial formula for arbitrary $T$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not possible here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo for arbitrary $T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The following function provides independent samples of the knock-out option at time T\n",
    "def sample_knockoutcall_value(T,a,b,S0,K,r,B,M):\n",
    "    value = S0*np.ones((M,1))\n",
    "    for t in np.arange(T-1,-1,-1):\n",
    "        random_values = bernoulli.rvs(size=(M,1),p=p_star)\n",
    "        up_moves = random_values*(1+b)\n",
    "        down_moves = -1*(random_values-1)*(1+a)\n",
    "        value = np.multiply(value,up_moves+down_moves)\n",
    "        value[value > B] = 0\n",
    "    return (np.maximum(value-K,0)/((1+r)**T))    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The last calculated value is 0.0000.\n",
      "The time needed for the last run was 8.30000 seconds.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEWCAYAAABfdFHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu4HFWd7vHvm01IAoRLSNSQADvRyEw8ct0iF2UQgaDjDDpyScgBBBxEccCHZ9BkELk4juIVeLxARBAYEFC8bDkyIVw9xwuwA4R7IMQwJEEIBEgEBkn4nT9q7aSz6e5d+1LdO9Xv53n66apVq2r9Vqc7a1fVqrUUEZiZmRVlWLMDMDOzcnNDY2ZmhXJDY2ZmhXJDY2ZmhXJDY2ZmhXJDY2ZmhXJDY9Ygkh6StH+z4zBrNPk5GrPmkPTjiPhEs+MwK5rPaMwaSNIOkr4jabO0/m5J32h2XGZF2qTZAZi1CklLgE8C1wNXAh3AK8BZafuPgZeBdmA/4GHgqIh4ovHRmg0en9GYNdda4I2K9RnAOcA2wCLgK80IymwwuaExa6ztgI8DRwO3AXOAWRXbfx4Rd0XEGuAqYNfGh2g2uHzpzKyxlkfEFQCSiIgHgNMrtv+5YvkVYItGBmdWBJ/RmDWJe5xZq3BDY2ZmhXJDY2ZmhfIDm2ZmViif0ZiZWaHc0JiZWaHc0JiZWaHc0JiZWaFa6oHNsWPHRnt7e7PDMDPbqMyfP/+5iBjX3/1bqqFpb2+nq6ur2WGYmW1UJD05kP2beulM0iGSFkpaJGlWle0jJF2btt8pqT2lt0t6VdJ96XVRo2M3M7N8mnZGI6kN+B5wELAUuFtSZ0Q8XJHtBOCFiHiHpOnAecCRadsTEeEBB83MhrhmntHsCSyKiMUR8VfgGuDQHnkOBS5Pyz8DPihJDYzRzMwGqJkNzQTgqYr1pSmtap40bPpLwLZp2yRJ90q6Q9L7axUi6URJXZK6VqxYMXjRm5lZLs1saKqdmfQcD6dWnqeBHSJiN+A04GpJW1YrJCLmRERHRHSMG9fvThNmZtZPzWxolgLbV6xPBJbXyiNpE2ArYGVEvBYRzwNExHzgCeCdhUdsZmZ91syG5m5giqRJkjYFpgOdPfJ0Asem5cOAWyMiJI1LnQmQNBmYAixuUNxmZtYHTet1FhFrJH0WmAu0AZdGxEOSzgW6IqIT+BFwpaRFwEqyxghgP+BcSWvI5lw/KSJWNr4WZmbWm5aaJqCjoyP8wKaZWd9Imh8RHf3d32OdmZlZodzQmJlZoXI1NJJ2lHRgWh4laXSxYZmZWVn02tBI+meyp/IvTkkTgV8WGZSZmZVHnjOak4F9gVUAEfE48JYigzIzs/LI09C8lsYiA9Y9ONk6XdXMzGxA8jQ0d0j6N2CUpIOAnwK/LjYsMzMrizwNzSxgBfAA8CngN8AXiwzKzMzKo9eRASLiDeCH6WVmZtYnvTY0kv5ElXsyETG5kIjMzKxU8ox1VjnswEjgcGBMMeGYmVnZ9HqPJiKer3gti4jzgQMaEJuZmZVAnktnu1esDiM7w/HIAGZmlkueS2ffqlheAywBjigkGjMzK508vc4+0IhAzMysnGo2NJJOq7djRHx78MMxM7OyqXdG4/swZmY2YDUbmog4p5GBmJlZOeXpdTYSOAF4F9lzNABExPEFxmVmZiWRZ6yzK4G3AdOAO8jmo1ldZFBmZlYeeRqad0TEmcDLEXE58PfAu4sNy8zMyiJPQ/N6en9R0v8CtgLaC4vIzMxKJc8Dm3MkbUM2NUAnsAVwZqFRmZlZadR7juatEfFMRFySkn4LeMRmMzPrk3qXzhZImifpeElbNSwiMzMrlXoNzQTgm8D7gcck/VLSkZJGNSY0MzMrg5oNTUSsjYi5EXEcsD1wGfBR4E+SrmpUgGZmtnHL0+uMiPgr8DDwCLAKmFpkUGZmVh51GxpJO0g6XdI9wA1AG3BoROzWkOjMzGyjV6/X2e/J7tP8FDgxIroaFpWZmZVGvedoZgO/jYhoVDBmZlY+9UZvvqORgZiZWTnl6gxgZmbWX25ozMysUJ7K2czMCpVnKuedgPeQDagJ8A9k456ZmZn1qtepnCXdBOweEavT+tlkXZ4HTNIhwAVkz+dcEhFf67F9BHAFsAfwPHBkRCxJ22aTzfy5FjglIuYORkw9/fLeZXxj7kKWv/gq2209itOn7cRHd5tQM70/+/hYPlYrHqvZ5bf6sRpJvfVelvQosEtEvJbWRwALIuJvBlSw1AY8BhwELAXuBmZExMMVeT4D7BwRJ0maDnwsIo6UNBX4CbAnsB1wM/DOiFhbr8yOjo7o6sr/ONAv713G7J8/wKuvrz/sqOFtfHyPCVw/f9mb0r/6T9l8cH3Zx8fysVrxWM0uv9WP1dfGRtL8iOjo006V++doaM4AjgB+AQTwMeC6iPiP/haajrs3cHZETEvrswEi4qsVeeamPH+QtAnwZ2AcMKsyb2W+emX2taHZ92u3suzFVzn3iTMZu/J/clQqvQ/Gk0c+lo9V5mM1u/wWO9ZzY0bypbd/GYAJW4/id7MO6FvxA2xoep34LCK+IulGslGcAY6LiHv7W2CFCcBTFetLgffWyhMRayS9BGyb0v/YY9+qTbSkE4ETAXbYYYc+Bbj8xVf7lH9QvlA+lo/VCsdqdvktfKw+/782CPLMsAmwGbAqIi6TNE7SpIj40wDLVpW0nh9nrTx59s0SI+YAcyA7o+lLgNttPYplL76a/SXw9vXpbRJrq5wJTtg6m0FhWZV/yFr7+Fg+Viseq9nlt/Kxttu68TO99PocjaSzgC+QDUkDMBz4z0EoeynZ9APdJgLLa+VJl862Albm3HfATp+2E6OGt22QNmp4GzPeu33V9NOn7dTnfXwsH6sVj9Xs8lv9WI2W54zmY8BuwD0AEbFc0uj6u+RyNzBF0iRgGTAdOKpHnk7gWOAPwGHArRERkjqBqyV9m6wzwBTgrkGIaQPdN8yq9dro2HFM3d4cfdnHx/KxWvFYzS7fx2qcPJ0B7oqIPSXdExG7S9oc+ENE7DzgwqUPA+eTdW++NN0POhfoiohOSSOBK8kaupXA9IhYnPY9AzgeWAN8LiJu7K28vnYGMDOzxvQ6+1eyM4aDgK+S/ef+k4i4sL+FNosbGjOzvmtEr7NvSjqIbGbNnYAvRcS8/hZoZmatpdeGRtJ5EfEFYF6VNDMzs7ryjN58UJW0Dw12IGZmVk71Rm/+NPAZ4O2S7q/YNBr4fdGBmZlZOdS7dHY1cCNZB4BZFemrI2JloVGZmVlp1Lx0FhEvpZGSLwBWRsSTEfEk8LqknkPFmJmZVZXnHs0PgL9UrL+c0szMzHqVp6FRVDxsExFvkH+MNDMza3F5GprFkk6RNDy9TgUWFx2YmZmVQ56G5iRgH7LxyLqH8j+xyKDMzKw88owM8CzZgJdmZmZ9lmeagHdKukXSg2l9Z0lfLD40MzMrgzyXzn5INhfN6wARcT8+wzEzs5zyNDSbRUTPuV7WFBGMmZmVT56G5jlJbydNlSzpMODpQqMyM7PSyPM8zMnAHOBvJC0D/gTMLDQqMzMrjTy9zhYDB6aZNYdFxOriwzIzs7LI0+tsW0kXAv8XuF3SBZK2LT40MzMrgzz3aK4BVgAfBw5Ly9cWGZSZmZVHnns0YyLiyxXr/y7po0UFZGZm5ZLnjOY2SdMlDUuvI4D/U3RgZmZWDnkamk+RTYL2WnpdA5wmabWkVUUGZ2ZmG788vc5GNyIQMzMrpzy9zk7osd4m6aziQjIzszLJc+nsg5J+I2m8pHcDfwR8lmNmZrnkuXR2lKQjgQeAV4AZEfG7wiMzM7NSyHPpbApwKnA9sAQ4WtJmBcdlZmYlkefS2a+BMyPiU8DfAY8DdxcalZmZlUaeBzb3jIhVABERwLckdRYblpmZlUXNMxpJnweIiFWSDu+x+bhCozIzs9Kod+mschbN2T22HVJALGZmVkL1GhrVWK62bmZmVlW9hiZqLFdbNzMzq6peZ4Bd0lhmAkZVjGsmYGThkZmZWSnUbGgioq2RgZiZWTnleY7GzMys35rS0EgaI2mepMfT+zY18h2b8jwu6diK9NslLZR0X3q9pXHRm5lZXzTrjGYWcEtETAFuSesbkDQGOAt4L7AncFaPBmlmROyaXs82ImgzM+u7ZjU0hwKXp+XLgWpTQ08D5kXEyoh4AZiHn98xM9voNKuheWtEPA2Q3qtd+poAPFWxvjSldbssXTY7U1LN53oknSipS1LXihUrBiN2MzPrgzxjnfWLpJuBt1XZdEbeQ1RJ635+Z2ZELJM0mmxU6aOBK6odJCLmAHMAOjo6/PyPmVmDFdbQRMSBtbZJekbS+Ih4WtJ4oNo9lqXA/hXrE4Hb07GXpffVkq4mu4dTtaExM7Pmatals06guxfZscCvquSZCxwsaZvUCeBgYK6kTSSNBZA0HPgI8GADYjYzs35QNvJ/gwuVtgWuA3YA/hs4PCJWSuoAToqIT6Z8xwP/lnb7SkRcJmlz4LfAcKANuBk4LSLW5ih3BfBkL9nGAs/1o1pl4Lq3rlaufyvXHfLVf8eIGNffAprS0AxlkroioqPZcTSD696adYfWrn8r1x0aU3+PDGBmZoVyQ2NmZoVyQ/Nmc5odQBO57q2rlevfynWHBtTf92jMzKxQPqMxM7NCuaExM7NCuaFJJB2Sph5YJOlNo0kPZZIulfSspAcr0qpOxaDMhame90vavWKfWtMy7CHpgbTPhd1jy+Wd7qFIkraXdJukRyQ9JOnUFqv/SEl3SVqQ6n9OSp8k6c4U27WSNk3pI9L6orS9veJYs1P6QknTKtKr/jZqldFoktok3SvphnpxlbTuS9J38z5JXSlt6H33I6LlX2QPfj4BTAY2BRYAU5sdVx/i3w/YHXiwIu3rwKy0PAs4Ly1/GLiRbCy5vYA7U/oYYHF63yYtb5O23QXsnfa5EfhQvTIaXPfxwO5peTTwGDC1heovYIu0PBy4M9XrOmB6Sr8I+HRa/gxwUVqeDlyblqem7/0IYFL6PbTV+23UKqMJn8FpwNXADfXiKmndlwBje6QNue9+wz+YofhKH+TcivXZwOxmx9XHOrSzYUOzEBiflscDC9PyxcCMnvmAGcDFFekXp7TxwKMV6evy1SqjyZ/Dr4CDWrH+wGbAPWRzOD0HbJLS132/yYZ22jstb5Lyqed3vjtfrd9G2qdqGQ2u80SyOa0OAG6oF1fZ6p7KXsKbG5oh992ve+lM0jBJR9TLUxK9TUmwMao1FUOtutZLX1olvV4ZTZEuhexG9ld9y9Q/XTq6j2xw2nlkf4W/GBFrUpbKmNfVM21/CdiWvn8u29Ypo5HOBz4PvJHW68VVtrpDNqL9TZLmSzoxpQ25737dhiYi3gA+Wy9PSdSbkqBsatW1r+lDiqQtyKaM+FxErKqXtUraRl3/iFgbEbuS/XW/J/C31bKl98Gqf9M/F0kfAZ6NiPmVyVWylq7uFfaNiN2BDwEnS9qvTt6m1TNPZ4B5kv5V2U3XMd2v/hQ2hC0Ftq9Ynwgsb1Isg+UZZVMwoA2nYqhV13rpE6uk1yujoZSN4n09cFVE/LyX2EpX/24R8SLZVBp7AVtL6p4GpDLmdfVM27cCVtL3z+W5OmU0yr7AP0paAlxDdvns/DpxlanuAETE8vT+LPALsj80htx3P09DczxwMtmIyfPTqyvHfhuTu4EpqSfJpmQ3CjubHNNA1ZqKoRM4JvVA2Qt4KZ36Vp2WIW1bLWmv1OPkmB7H6m26h0KlmH4EPBIR367Y1Cr1Hydp67Q8CjgQeAS4DTisSmyVMR8G3BrZhfZOYHrqmTUJmEJ2I7jqbyPtU6uMhoiI2RExMSLaU1y3RsTMOnGVpu4AkjZXNvkjyka1P5hsypSh991vxg2sofgi65HxGNn17TOaHU8fY/8J8DTwOtlfISeQXUe+BXg8vY9JeQV8L9XzAaCj4jjHA4vS67iK9I70BX4C+C7rR5SoWkaD6/4+stP5+4H70uvDLVT/nYF7U/0fBL6U0ieT/We5CPgpMCKlj0zri9L2yRXHOiPVcSGpd1G930atMpr0G9if9b3OWqLuKYYF6fVQd3xD8bvf6xA06bLEp8m60EJ2an5xRLxed0czMzNyjHUm6RKy/vmXp6SjgbWRJiczMzOrJ09DsyAiduktzczMrJo8nQHWSnp794qkyUCv0yabmZlB9nRsb04HbpO0mOxm0o7AcYVGVZCxY8dGe3t7s8MwM9uozJ8//7mIGNff/es2NJKGAa+SdffbiayheTQiXutvgT2OfwhwAdmYQpdExNd6bB8BXAHsATwPHBkRS9IT4I+Q9RAB+GNEnNRbee3t7XR1la1ntplZsSQ9OZD96zY0EfGGpG9FxN5k3ScHjaQ2sq52B5F1yb1bUmdEPFyR7QTghYh4h6TpwHnAkWnbE5E9DW1mZkNYnns0N0n6ePfw0INoT2BRRCyOiL+SPdl7aI88h7K+t9vPgA8WEIeZmRUoT0NzGtkDSa9JWiVptaR6Y0nllWcgy1qD4AFMUjYHxR2S3l+rEEknSuqS1LVixYpBCNvMzPqit9GbBbwrIoZFxKYRsWVEjI6ILQeh7DwDttXK8zSwQ0TsRpqLQlLVmCJiTkR0RETHuHH9vpdlZmb91NvozUE2UFsR8gxkWXUQvIh4LSKeTzHOJxse4Z0FxWlmZgOQ59LZHyW9p4Cy8wxkWXUQvDSQYBuse65nCtmscGZmNsTkeY7mA8CnUve2l8kuZ0VE7DyQgiNijaTPko0c2gZcGhEPSToX6IqITrJRea+UtIhsOO/paff9gHMlrSF7ePSkiFg5kHjMzKwYeYag2bFaekQMqF91M3R0dISfozEz6xtJ8yOio7/717x0JukAWNegDIuIJ7tfZA9QmpmZ9arePZpvVixf32PbFwuIxczMSqheQ6May9XWzczMqqrX0ESN5WrrZmZmVdXrdTZZUifZ2Uv3Mml9UuGRmZlZKdRraCrHHftmj209183MzKqq2dBExB2NDMTMzMopz8gAZmZm/eaGxszMCpW7oZG0eZGBmJlZOfXa0EjaR9LDZFMnI2kXSd8vPDIzMyuFPGc03wGmAd3D8i8gG9TSzMysV7kunUXEUz2S1hYQi5mZlVCeaQKekrQPEGnemFNIl9HMzMx6k+eM5iTgZGAC2YyXu6Z1MzOzXtU9o0mzWB4dETMbFI+ZmZVM3TOaiFjLhkPRmJmZ9UmeezS/k/Rd4FqyqZwBiIh7CovKzMxKI09Ds096P7ciLYADBj8cMzMrm14bmoj4QCMCMTOzcspzRoOkvwfeBYzsTouIc2vvYWZmlskzBM1FwJHAv5BNenY4sGPBcZmZWUnkeY5mn4g4BnghIs4B9ga2LzYsMzMrizwNzavp/RVJ2wGv46mczcwspzz3aG6QtDXwDeAesh5nlxQalZmZlUaeXmdfTovXS7oBGBkRLxUblpmZlUWvDY2kY6qkERFXFBOSmZmVSZ5LZ++pWB4JfJDsEpobGjMz61WeS2f/UrkuaSvgysIiMjOzUsk18VkPrwBTBjsQMzMrpzz3aH5N1tMMsoZpKnBdkUGZmVl55LlH882K5TXAkxGxtKB4zMysZPLco7mjEYGYmVk55bl0tpr1l8422ARERGw56FGZmVlp5Ll09h3gz2Q9zQTMBEZHxNeLDMzMzMohT6+zaRHx/YhYHRGrIuIHwMeLDszMzMohT0OzVtJMSW2ShkmaCawdjMIlHSJpoaRFkmZV2T5C0rVp+52S2iu2zU7pCyVNG4x4zMxs8OVpaI4CjgCeAZ4lm4/mqIEWLKkN+B7wIbIu0zMkTe2R7QSy6QneQXYJ77y071RgOtlkbIcA30/HMzOzISZPr7MlwKEFlL0nsCgiFgNIuiaV83BFnkOBs9Pyz4DvSlJKvyYiXgP+JGlROt4fCoiTc379EA8vX1XEoc3MCjd1uy056x/e1bTya57RSPpnSVPSsiRdKuklSfdL2n0Qyp4APFWxvjSlVc0TEWuAl4Btc+7bXY8TJXVJ6lqxYsUghG1mZn1R74zmVODHaXkGsAswGdgNuAB4/wDLVpW0nt2oa+XJs2+WGDEHmAPQ0dFRNU9vmvmXgJnZxq7ePZo1EfF6Wv4IcEVEPB8RNwObD0LZS9lwSuiJwPJaeSRtAmwFrMy5r5mZDQH1Gpo3JI2X1D01wM0V20YNQtl3A1MkTZK0KdnN/c4eeTqBY9PyYcCtEREpfXrqlTaJbJDPuwYhJjMzG2T1Lp19CegC2oDOiHgIQNLfAYsHWnBErJH0WWBuKuPSiHhI0rlAV0R0Aj8Crkw3+1eSNUakfNeRdRxYA5wcEYPS5drMzAaXshOEGhuzy1WjI+KFirTN035/aUB8g6qjoyO6urqaHYaZ2UZF0vyI6Ojv/nW7N6eeXi/0SHu5v4WZmVnr6c/EZ2ZmZrm5oTEzs0LlGb0ZSROAHSvzR8RviwrKzMzKI898NOcBR5L18Oru2RWAGxozM+tVnjOajwI7pXHFzMzM+iTPPZrFwPCiAzEzs3LKc0bzCnCfpFuAdWc1EXFKYVGZmVlp5GloOnnz0DBmZma55JmP5vJGBGJmZuWUp9fZFOCrZLNgjuxOj4jJBcZlZmYlkaczwGXAD8gGr/wAcAVwZZFBmZlZeeRpaEZFxC1kA2k+GRFnAwcUG5aZmZVFns4A/yNpGPB4GtZ/GfCWYsMyM7OyyHNG8zlgM+AUYA/gf7N+MjIzM7O68vQ6uxtAUkTEccWHZGZmZdLrGY2kvSU9DDyS1neR9P3CIzMzs1LIc+nsfGAa8DxARCwA9isyKDMzK49c89FExFM9ktZWzWhmZtZDnl5nT0naBwhJm5J1Cnik2LDMzKws8pzRnAScDEwAlgK7pnUzM7Ne5el19hwwswGxmJlZCdVsaCRdWG9HTxNgZmZ51DujOQl4ELgOWA6oIRGZmVmp1GtoxgOHA0eSDah5LXB9RLzQiMDMzKwcanYGiIjnI+KiiPgA8Alga+AhSUc3KjgzM9v45ZmPZndgBnAQcCMwv+igzMysPOp1BjgH+AjZMzPXALMjYk2jAjMzs3Kod0ZzJrAY2CW9/kMSZJ0CIiJ2Lj48MzPb2NVraCY1LAozMyutmg1NRDzZyEDMzKyccg2qaWZm1l9uaMzMrFB5Jj47NU+amZlZNXnOaI6tkvaJQY7DzMxKqt5zNDOAo4BJkjorNo0mzbZpZmbWm3rdm38PPA2MBb5Vkb4auH8ghUoaQzZ2WjuwBDii2hhqko4FvphW/z0iLk/pt5ONxfZq2nZwRDw7kJjMzKwY9cY6ezIibo+IvYFHyc5kRgNLB2GEgFnALRExBbglrW8gNUZnAe8F9gTOkrRNRZaZEbFrermRMTMbovJ0BjgcuItsJOcjgDslHTbAcg8FLk/LlwMfrZJnGjAvIlams515wCEDLNfMzBqs10E1yS5dvaf7rEHSOOBm4GcDKPetEfE0QEQ8LektVfJMAJ6qWF+a0rpdJmktcD3ZZbWoVpCkE4ETAXbYYYcBhGxmZv2Rp6EZ1uPS1PPkOxO6GXhblU1n5Iyt2kRr3Y3JzIhYJmk0WUNzNHBFtYNExBxgDkBHR0fVxsjMzIqTp6H5L0lzgZ+k9SOB3/S2U0QcWGubpGckjU9nM+OBavdYlgL7V6xPBG5Px16W3ldLuprsHk7VhsbMzJqr1zOTiDgduBjYmWwU5zkR8YUBltvJ+udzjgV+VSXPXOBgSdukTgAHA3MlbSJpLICk4WRTGTw4wHjMzKwgec5oAH4HvE526equQSj3a8B1kk4A/pusowGSOoCTIuKTEbFS0peBu9M+56a0zckanOFAG9n9oh8OQkxmZlYA1biHvj6DdATwDbLLVgLeD5weEQPpDNAUklYAvY1KPRZ4rgHhDEWue+tq5fq3ct0hX/13jIhx/S0gT0OzADioZ6+ziNilv4UOZZK6IqKj2XE0g+vemnWH1q5/K9cdGlP/PGOd9avXmZmZGfS/19mNxYVkZmZl0mtDExGnS/on4H1k92jmRMQvCo+seeY0O4Amct1bVyvXv5XrDg2of6/3aN60g9QGTI+Iq4oJyczMyqTmvRZJW0qaLem7kg5W5rPAYrIxz8zMzHpV84xG0q+AF4A/AB8EtgE2BU6NiPsaFqGZmW3cIqLqC3igYrmNrNEZXSv/xv4iGxl6IbAImNXsePoY+6Vkw/g8WJE2hmzE68fT+zYpXcCFqZ73A7tX7HNsyv84cGxF+h7AA2mfC1n/B0rVMhpc9+2B24BHgIfI/hBqpfqPJHuIekGq/zkpfRJwZ4rtWmDTlD4irS9K29srjjU7pS8EpvX226hVRhM+gzbgXuCGFqz7kvTdvA/oGqrf/XoVuKfeeple6Yv6BDCZ7KxtATC12XH1If79gN3ZsKH5evcPg2y+n/PS8ofJeg0K2Au4s+KLszi9b5OWu7+gdwF7p31uBD5Ur4wG13189w+GbL6kx4CpLVR/AVuk5eHpP7+9gOvI7qUCXAR8Oi1/BrgoLU8Hrk3LU9P3fgTZf6JPpN9Fzd9GrTKa8BmcBlzN+oamleq+BBjbI23IfffrVWAtsCq9VgNrKpZXNeNDLfAfa29gbsX6bGB2s+PqYx3a2bChWQiMT8vjgYVp+WJgRs98wAzg4or0i1PaeODRivR1+WqV0eTP4VfAQa1Yf2Az4B6yyQKfAzZJ6eu+32RjCO6dljdJ+dTzO9+dr9ZvI+1TtYwG13ki2eSJBwA31IurbHVPZS/hzQ3NkPvu15thsy0itkyv0RGxScXylrX220j1NvfNxmiDOX+A7jl/atW1XvrSKun1ymgKSe3AbmR/1bdM/SW1SbqP7PLpPLK/wl+M9TPhVsa8rp5p+0vAtvT9c9m2ThmNdD7weeCNtF4vrrLVHbLxJ2+SND/NvQVD8Lufd1DNsqs3903Z1KprX9OHFElbkM1N9LmIWCVVCzvLWiVto65/RKwFdpW0NfAL4G+rZUvvfa1ntT9Gh8TnIukjwLMRMV/S/t3JVbKWru4V9o2I5WnyyHmSHq2Tt2nffQ8lk1lKdlO520RgeZNiGSzPpLl+6DH88bX8AAACp0lEQVTnT6261kufWCW9XhkNlUbyvh64KiJ+3ktspat/t4h4kWzw272ArSV1/yFZGfO6eqbtWwEr6fvn8lydMhplX+AfJS0BriG7fHZ+nbjKVHcAImJ5en+W7I+MPRmC3303NJm7gSmSJknalOxGYWeTYxqoWnP+dALHpOei9gJeSqe+Vef/SdtWS9pL2WnCMT2OVa2Mhkkx/Qh4JCK+XbGpVeo/Lp3JIGkUcCBZD7zbgMOqxFYZ82HArZFdaO8EpksaIWkSMIXsRnDV30bap1YZDRERsyNiYkS0p7hujYiZdeIqTd0BJG2eZhkmTZ9yMNncXEPvu9+MG1hD8UXWI+MxsuvbZzQ7nj7G/hPgabI5g5YCJ5BdR76FrPvhLcCYlFfA91I9HwA6Ko5zPFk3xkXAcRXpHekL/ATwXdZ3caxaRoPr/j6y0/n7ybp43pf+LVul/juTde29P8X4pZQ+mew/y0XAT4ERKX1kWl+Utk+uONYZqY4LSb2L6v02apXRpN/A/qzvddYSdU8xLGB91/Yz6n0vm/nd7/MQNGZmZn3hS2dmZlYoNzRmZlYoNzRmZlYoNzRmZlYoNzRmZlYoNzRmDSTp9+ld6f3synWzMnL3ZrMmkHQa2SC1OwF/Be6IiJuaG5VZMXxGY9ZAkv4CENkoBmOBU4D/ioibJO0v6XZJP5P0qKSrfKZjZeBBNc2aQNLnyMbMuhA4RNJIspEddgPeRTam1O/IxvP6f82K02ww+IzGrDkuiIhLgJcj4gzg5pR+V0QsjYg3yIbTaW9WgGaDxQ2NWRNEujkaEWdXrgOvVWRbi686WAm4oTEzs0K5oTEzs0K5e7OZmRXKZzRmZlYoNzRmZlYoNzRmZlYoNzRmZlYoNzRmZlYoNzRmZlYoNzRmZlao/w+rR+nZFqCWcwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Value of pi via Monte Carlo\n",
    "n=10000\n",
    "mc_values = 50\n",
    "tic_2 = time.time()\n",
    "x_values = []\n",
    "pi_M = []\n",
    "pi_M_RMSE = []\n",
    "for i in range(1,mc_values,1):\n",
    "    M=n*i\n",
    "    if i == mc_values-1:\n",
    "        tic_2 = time.time()    \n",
    "    x_values.append([M])\n",
    "    values_call = sample_knockoutcall_value(T,a,b,S0,K,r,B,M)\n",
    "    pi_M.append([np.mean(values_call)])\n",
    "    pi_M_RMSE.append([np.sqrt(np.var(values_call,ddof=1))/np.sqrt(n*i)])\n",
    "toc_2 = time.time()\n",
    "print( \"The last calculated value is %.4f.\" %(float(pi_M[len(pi_M)-1][0])))\n",
    "print(\"The time needed for the last run was %.5f seconds.\" %(toc_2-tic_2))  \n",
    "#The following will plot the calculated value for pi\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(x_values, pi_M, 'o')\n",
    "plt.plot(x_values, [pi_1]*len(x_values), '-')\n",
    "plt.plot(x_values, pi_1-1.96*np.array(pi_M_RMSE), '-')\n",
    "plt.plot(x_values, pi_1+1.96*np.array(pi_M_RMSE), '-')\n",
    "plt.title('i*n')\n",
    "plt.ylabel('Expected Value')\n",
    "#The following will plot the root mean square error\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(x_values, pi_M_RMSE, '-')\n",
    "plt.xlabel('i*n')\n",
    "plt.ylabel('Root Mean Square Error')\n",
    "plt.show()    "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}