FF3F.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0.5\n"
     ]
    }
   ],
   "source": [
    "# Pandas to read csv file and other things\n",
    "import pandas as pd\n",
    "# Datareader to download price data from Yahoo Finance\n",
    "import pandas_datareader as web\n",
    "# Statsmodels to run our multiple regression model\n",
    "import statsmodels.api as sm\n",
    "# To download the Fama French data from the web\n",
    "import urllib.request\n",
    "# To unzip the ZipFile \n",
    "import zipfile\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize,root_scalar\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "print(pd.__version__)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_fama_french():\n",
    "    # Web url\n",
    "    ff_url = \"https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_CSV.zip\"\n",
    "    # Download the file and save it\n",
    "    \n",
    "    # We will name it fama_french.zip file\n",
    "    urllib.request.urlretrieve(ff_url,'fama_french.zip')\n",
    "    zip_file = zipfile.ZipFile('fama_french.zip', 'r')\n",
    "    \n",
    "    # Next we extact the file data\n",
    "    zip_file.extractall()\n",
    "    \n",
    "    # Make sure you close the file after extraction\n",
    "    zip_file.close()\n",
    "    \n",
    "    # Now open the CSV file\n",
    "    ff_factors = pd.read_csv('F-F_Research_Data_Factors.csv', skiprows = 3, index_col = 0)\n",
    "\n",
    "    # We want to find out the row with NULL value\n",
    "    # We will skip these rows\n",
    "    # need to convert series to array to call \"nonzero\" function\n",
    "    ff_row = np.array(ff_factors.isnull().any(1)).nonzero()[0][0]\n",
    "    \n",
    "    # Read the csv file again with skipped rows\n",
    "    ff_factors = pd.read_csv('F-F_Research_Data_Factors.csv', skiprows = 3, nrows = ff_row, index_col = 0)\n",
    "    \n",
    "    # Format the date index\n",
    "    ff_factors.index = pd.to_datetime(ff_factors.index, format= '%Y%m')\n",
    "    \n",
    "    # Format dates to end of month\n",
    "    ff_factors.index = ff_factors.index + pd.offsets.MonthEnd()\n",
    "    \n",
    "    # Convert from percent to decimal\n",
    "    ff_factors = ff_factors.apply(lambda x: x/ 100)\n",
    "    return ff_factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            Mkt-RF     SMB     HML      RF\n",
      "2020-05-31  0.0558  0.0247 -0.0495  0.0001\n",
      "2020-06-30  0.0246  0.0270 -0.0222  0.0001\n",
      "2020-07-31  0.0577 -0.0218 -0.0131  0.0001\n",
      "2020-08-31  0.0763 -0.0026 -0.0295  0.0001\n",
      "2020-09-30 -0.0363  0.0010 -0.0259  0.0001\n"
     ]
    }
   ],
   "source": [
    "ff_data = get_fama_french()\n",
    "print(ff_data.tail())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Last day of FF data\n",
    "ff_last = ff_data.index[ff_data.shape[0] - 1].date()\n",
    "# Build the get_price function\n",
    "# We need 3 arguments, ticker, start and end date\n",
    "def get_price_data(ticker, start, end):\n",
    "    price = web.get_data_yahoo(ticker, start, end)\n",
    "    price = price['Adj Close'] # keep only the Adj Price col\n",
    "    return price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Date\n",
      "2020-09-24    16.030001\n",
      "2020-09-25    16.330000\n",
      "2020-09-28    16.549999\n",
      "2020-09-29    16.570000\n",
      "2020-09-30    16.680000\n",
      "Name: Adj Close, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "# Get Price data for Fidelity's fund\n",
    "price_data = get_price_data(\"FCNTX\", \"1980-01-01\", \"2020-09-30\")\n",
    "# Make sure to only have data upto last date of Fama French data\n",
    "price_data = price_data.loc[:ff_last]\n",
    "print(price_data.tail())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            portfolio\n",
      "Date                 \n",
      "2020-05-31   0.064588\n",
      "2020-06-30   0.040446\n",
      "2020-07-31   0.071716\n",
      "2020-08-31   0.096310\n",
      "2020-09-30  -0.048488\n"
     ]
    }
   ],
   "source": [
    "def get_return_data(price_data, period = \"M\"):\n",
    "    \n",
    "    # Resample the data to monthly price\n",
    "    price = price_data.resample(period).last()\n",
    "    \n",
    "    # Calculate the percent change\n",
    "    ret_data = price.pct_change()[1:]\n",
    "    \n",
    "    # convert from series to DataFrame\n",
    "    ret_data = pd.DataFrame(ret_data)\n",
    "    \n",
    "    # Rename the Column\n",
    "    ret_data.columns = ['portfolio']\n",
    "    return ret_data\n",
    "    \n",
    "ret_data = get_return_data(price_data, \"M\")\n",
    "print(ret_data.tail())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            portfolio  mkt_excess     SMB     HML      RF  port_excess\n",
      "2020-05-31   0.064588      0.0558  0.0247 -0.0495  0.0001     0.064488\n",
      "2020-06-30   0.040446      0.0246  0.0270 -0.0222  0.0001     0.040346\n",
      "2020-07-31   0.071716      0.0577 -0.0218 -0.0131  0.0001     0.071616\n",
      "2020-08-31   0.096310      0.0763 -0.0026 -0.0295  0.0001     0.096210\n",
      "2020-09-30  -0.048488     -0.0363  0.0010 -0.0259  0.0001    -0.048588\n"
     ]
    }
   ],
   "source": [
    "# Merging the data\n",
    "all_data = pd.merge(pd.DataFrame(ret_data),ff_data, how = 'inner', left_index= True, right_index= True)\n",
    "# Rename the columns\n",
    "all_data.rename(columns={\"Mkt-RF\":\"mkt_excess\"}, inplace=True)\n",
    "# Calculate the excess returns\n",
    "all_data['port_excess'] = all_data['portfolio'] - all_data['RF']\n",
    "print(all_data.tail())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9eZQc93Xf+6mt99kHgxlgMAPMglUkARAgoH0zSVlyFC2OIlsv1rEix7KdyHGip6ecnOPEju3DvJNEii3Jsf0U2XpPlmXHjmjH1C45pi0SIEASJBYCs+/79L7V9nt/9FSze6Z7prqnZzhD9vccHAAzVb+urq7+1q17v/d7JSEEddRRRx117A7kV/oA6qijjjpeS6iTbh111FHHLqJOunXUUUcdu4g66dZRRx117CLqpFtHHXXUsYtQt/h9XdpQRx111FE5pHK/qEe6ddRRRx27iDrp1lFHHXXsIuqkW0cdddSxi6iTbh111FHHLqJOunXUUUcdu4g66dZRRx117CLqpFtHHXXUsYuok24dddRRxy6iTrp11FFHHbuIOunWUUcddewi6qRbRx111LGLqJNuHXXUUccuok66ddRRRx27iK1cxuqooyyEENi2TTabxTRNVFVFlmUURUGWZWRZRpLKmi3VUcdrEtIWgynr1o51bIAQAsuyME2z6N8AqVSK2dlZBgYGAPIk7Pypk3EdrxGUvcDrkW4drrGebCVJQpZlbNvO/1tRFHRdR1EUhBAIITAMA13Xi4i2TsZ1vFZRJ906toQQAtM0sSyriGxLoZA0JUkqSaLO05VpmhiGUfS7OhnX8WpHnXTrKAshBLquMzs7S2dn56Zk60CSJLZIWeUJdD2RrifjqakpmpqaaGxsRJIkFEXJ540dcq6TcR37DXXSrWMDbNvOpxFs22ZiYoJDhw6V3X59dLsV6W61jvO3k8ZwiH59asP5XbnIuE7IdexF1Em3jjxs286nEYB8dFkJiW6HdDdbs/BvB87rFBbyCveRZRlVVetkXMeeQp10X+MoLHbZtg2Uz8W6QS1Jd6u1tiJjR84GMDw8zLFjx4pkbYWpijoZ17FbqJPuaxSOxtZJIcD2yNbBTkS61RxD4d8AyWQyT65O+kTX9aL9CtMUTnRcJ+M6ao066b7G4CgRYrEYgUAgTyq1Ipa9QLqbYTNFhUPGTs7YQamccV1RUUe1qJPuawSFGlvDMLh16xYPPfRQzYljN9MLlcCNosItGU9MTNDd3Y2maXV5Wx0Vo066r3KUamhQVTXf0FBr7OVIt5r3W4qMw+EwPT099caPOqpCnXRfpdiqoWGniPG1QCybNYjUGz/q2Ap10n2VwSFbR0JVihx2IqWwE9iPUbPbxg8Hq6urhEIhAoFAvfHjNYI66b5KUNjQALVRIryaUEvyXl9oc4NyZBwOh1FVFa/XW2/8eI2gTrr7HI7sa2lpiUgkQn9//65/IR2t715HLc9LrdYSQuQJdf3Pod748WpEnXT3Ico1NBiGsWNfulLRnRCC+fl5xsbG8oW5QCBANptlaWmJYDCI3+/fE40WexXlCpqVNH44sCyLTCZDS0tLvfFjD6NOuvsImzU0VNquWwkc8nO+uLZtMzc3x8TEBC0tLZw9ezYfqaVSKaLRKPF4nPn5edLpdJ6Mg8EgoVCIYDCIz+erEwGVpyrKkTFAOp1mamqKYDBYb/zYw6iT7j7AerIt1dAgy3LeM6HWcEjXtm1mZmaYnJykvb2dBx98EK/Xi23b6LqOLMuEQiE0TaOvry+/v23bpFIpEokE0WiU2dlZMpkMsiznydj581ojY9u2t3Ruq2Qth1QLUUpr7Pxdb/zYfdRJdw+jlMa2XHTimInvBCRJYnJyktnZWTo6Orh48SIej8f1/g4Zh0Khop9blkUqlSKZTJYkY8MwME2TxsZGvF7vtohgr6YpqinKlUM5Aq+mC0+WZYQQaJqGpml1Mq4h6qS7B+GQ7dTUFO3t7fmLfjNUQ7pbfeFN02RqaopYLEZLSwuXLl1CVUtfMtV8GRVFoaGhgYaGhqKfO2Q8MTFBKpXi7t27ZLNZZFkuioqDwWBFZLwXCaPWkW4la21GxkIIhoaGaGtro6WlJf+7utZ4+6iT7h7C+oaG+fl5WlpaXEWVlZKuE8mU+rKYpsnExARzc3McPnyYlpYWent7yxJureGQcWNjIx6Ph87OTiBHxslkkmQySTgcZnp6mmw2i6IoBAKBfL44GAzi8Xj2BRHsRqRbKRwyFkLg8Xjy6Yp640dtUCfdPYByDQ2KorjO0yqKUhHpOm5bhV9SXdeZmJhgcXGR7u5uXv/616MoCuFweMu1a0ke5aAoCo2NjTQ2Nhb93DTNfM54ZWWFycnJ/Jw2h4RN0ySbzW47TVFrvJKR7lawLKsoP1xp4weQJ99S8rbXKuqk+wpiq4aGSqLXSgtphWtns1nGx8dZXl6mp6eH17/+9UVfXjfSrVfyS6SqalkydiJj0zS5c+cOhmGgquqGNEUlOepaYi9Gug4sy3K13lZk7Fzjd+7cYXBwME++6z2NXyuKijrpvgIoNaGh1IVWSaRbTXohnU4zMjJCOBymt7eXwcHBkl+yrdbeq23AqqrS1NREU1MTMzMznDt3DgDDMPJkvLS0xPj4eEkydpQYO4m9TLqOGqJarCfjTCaDpmn5/79WGz/qpLtLqGZCQyWkW8kFmU6nSSQS3Lx5k76+Pk6ePLnp/q+GJoXC49c0jebmZpqbm4u2WU/GY2NjmKaJpmlFZPxKtxSXw06nF7aLwnl3UFnjR2HKbb83ftRJd4dRqLGdm5sjlUrR19fn6kKptQwsmUwyOjpKMpnE5/Nx5syZDTKuWh3HbuR4a41yZKzrep6MFxYWSKfTXL16tYiMnSJeNcXGWpJuLYudtSbxrbBZ40ehvM1p/HAkjF6vN1/Ac4Ka9ammvYQ66e4QSjU0OJGr2y9ZJZHuZkgkEoyMjJDJZOjr66O9vZ0XX3zRdcRWaaS7voOtWtQywt7OsXg8HjweT146FYvFuHjxYhEZz8/P53PHHo9nQ854N5QftSbJvXTjLBXRLi0t0dnZmSdjgG9/+9vcuHGD3/qt33olDtMV6qRbY2zW0FApiW430o3FYoyMjGCaJn19fbS2tuYv3ErWrpZ0X+1YT8aQ+/wLyXhubo5kMollWXg8niJZWyAQqOnx1Jp0a0m4O3E9WJaFqqpFxxmNRmlqaqr5a9USddKtEUqR7fovgKqqFZFutZFuJBJhZGQEgL6+viJScFCpMqKa9MJrEZIk4fV68Xq9tLa25n++noxnZmby3Xg3btzYEBlXk0vdiUi3VtipVMX6G0M0Gt2QHtprqJPuNrHVhIZCVKMwWG9cshlM0+TatWvIsszAwMCmd/ydjnRrgb0YMVd7POXI+OrVq5w4caKIjJPJJLZt4/P5iog4EAhsSsa7nYOtBLUuypVDNBqlt7d3x19nO6iTbpUQQpBOp4lEIrS0tGxKtg4qjVzdNDwIIVhdXWVkZARd1zl79qyrIoLTkeYGld4s9iJZ7lVIkoTP58Pn89HW1pb/uRCCTCaTJ+PV1VVSqdSmZFxL0q11PncnlBClUI90X4UobGhIpVJMT08XfVk2g6qqG3SJm2EzkhZC5GVNfr+f06dPc/PmTYLBoKu1nY40t9vWc7o57FZxSZIk/H4/fr+f9vb2otcvR8ZON1g6nc6TcbUkXOuo2TTNmhYTy5F4LBark+6rBaUaGjRN29HCWKnthRAsLCwwNjZGQ0MD9913X74g40Q7biKKSnO65Ug0lUoxOjrK6upqXkKVyWQIh8OufSPKoZIbw27hla7ob0bGN27cIBQKkU6nWV5eJpVKIYTA7/dviIy3ItRaR6a1Xq8cicdisXohbT9jq4YGRVFqFrlutb0Qgrm5OcbHx2lububs2bP4/f6i7Z1WYDddVJXmdNdv65BtIpGgr6+PgYEBbNsmmUwSi8VYXV1ldnYWwzDyZFxYud8t85xao5akW8unASe9deDAAbxeb9FrpNPpfGTskDGwgYz9fn+ejPcD6ZZaLxqNliwc7yXszyt/h7HZhIZCVEqilX5ZHRKdnp5mcnKS1tZWzp8/j8/nK7v9TigSCtMF68n2zJkzQK6BwGm79fv9HD16NH9TcKr2iUSiSEJVmJt0JuLu1UKQg73ajQalUwKSlJvaEQgEOHDgQNG2hWS8uLhIOp0GcmTs8XjQdZ1UKrWtkUsOdoLES9246zndfQZH9mVZVtkJDYWopBhVKWzbZnFxkcXFRXw+HxcuXNjyUX2nFAmyLJPJZLh582YR2TrnZf0669cup2fNZrMkEgmSySQrKysbIjCHyGtBTrWOKmuB3SDdcij0Jl6/RjqdZmlpCdu2GR0dzZPx+ikflZBxrXO65dZz+6T3SqJOurxMtpFIhEgkQnd39yvW0+2Yl8/MzNDW1kZTUxPHjx93te9OuJKlUinm5ubIZDKcOnWqiGzLwa0rmVO1L8xNFkZg8XicaDTKyspKfpqEk6IIhUKviGfuXvZKKDfkshI4ZKzrOoZhMDg4mF/b0RVXM/+uXGRaLUqR7n4p3r6mSXd9Q4NlWUSjUY4cObLrx+JMaZiZmeHQoUNcunQJ27a5ceOG6zUq8dSVZXmD92khCtMITU1NdHR00NHR4Wrt7agXCiMw57M5cuRIfppEIpHYYGBemCsu5wy2F6PTWke66w1ltoP1N4RyI5fczr/LZDI19UPYjHT3SutyObwmSbecabimaRUVxrbz+s6FYRgGk5OTzM3NFRmHA/lUh1tU4qlbLip27B4L0whzc3MbXJ/Wo/BCr5VkrHCdcqN9HM/cRCJR5Azm+B84hFyrKGgvR7q1hNscrNv5d8vLyywvLzM1NbVhykc1xvKmaW6obWQymZq3Vu8EXlOkW4psCz/sSnW0zhqVfHmc4ptt2/kpDUeOHCkiWwe1kJi53bYU2Rb6NOxVnW6hZ64DR3Hi5IudlturV6/i8/k2+B9UQnx7OdKtJbZb+Fp/kzRNk87OTkKhUMknlkrn35WKdCORyJ52F3PwmiDdrSY0OKjG68Ahard6VEmSuHv3LpFIhN7e3g1TGtZvWwmqKaRtRraF21bTkVaLinU15C1JEh6Ph9bW1nzLbSKR4MKFC/nGgkQiUSSfKnwUDoVCZUfB7wc/glpgp7x5NxtGWm7+XakpH6VIdz8oF+BVTrpuJzQ4qIZ03e6TyWQYGxsjFotx4MABTp06VfMvXCWkaxgGi4uLrK6uliVbB24i18KobSUe5je+8QWG5sZpa2zh1//Jp3hw8L6K3ovzurVEucYCp3iXSCSIx+P5wqETfRVGxq+lSLeWI4y2ugFXOv8um81iGAaxWCyv+w6Hw3u+MQJehaRb2NDgdG21tbW5urir+QJs5RyWSqXyZHvs2DF0Xae9vX1HIhw3pOtEtpFIBL/fz/nz57d835WQuRCC//sv/4D56BIdLe0kMyn+9R/8B77+b77AwZYDWy/wCqCcfKow+nK+8JlMBsMwuHfvXlFkXE1lfq+T7m50kG2FcvPvnnvuOTo7O9F1nZWVFb74xS/y/e9/H9u2+cVf/EXOnDnDhz70oSJt8np861vf4pd/+ZexLIuPf/zjfOYznyn6/UsvvcTP/uzP8uyzz/Kbv/mbfOpTn3K976bvyfWWexylGhqcfvSdvLDLdaUVTmno6+vj9OnTSJLEyspKTYzJS2GzQtr6NEJ3dzfT09Oub0aFka5lWUwuzaLICt3tnUV63UQmxdTyLAdb2pEkiZA/SDgeYWh2rGLSfaU9HEpFX4lEgvHxcTo6OkgkEiwsLDAyMoJlWXi93g3NHrvlClbr87TXR//Ytk1LS0t+zc9+9rN8/etfZ3JykkcffZRbt25tWp+xLItf+qVf4rvf/S7d3d1cvHiR9773vZw+fTq/TWtrK7/927/NN77xjYr33Qz7nnQ3a2jYDTXC+kg3Ho8zOjpKJpOhv79/Q5Rdq2kQpeCMKylEuZxtPB6vynshnkrwi5//t9yauAdC8IbTF/gvP/+reYL0e7xoqopuGmgeTy6fbts0B/f+Y58bCCFQFGXDWB/HL9cp3k1NTeWNaJxmDydN4bTb7mUlxE7PR9suSr3fWCxGV1cXb3zjG3njG9+46f5Xr15lYGCAvr4+AD784Q/z+OOPFxGnI5P867/+64r33Qz7lnQ3m9DgQFXVLaVO61GtGqFwSkN/f3+RZ2qp7XfieArTAFsVyCqJIgsLab/z+Jd5YfQO7U259/fkzat89Qf/k4+8/f0AqIrKLzzyEX73u39M2sgiELzn4js40+uuwWOvo9w5K/TLXW/R6DR7JBKJonZbp9tueXmZUChUlXTKwV4n3Z1AKQPzrq4uV/vOzMwU6fG7u7u5cuXKju8L+5B03UxocKCqKslksqL1K1Uj6LrO0NAQPp+v7JSGQlQ7ssct6Waz2bLtuqXWdYNCgr49eQ+/9+XKvqqq3JkcLtr+zacf4kzfSVbSUdoaW3hw4L6qyOSVTi+UQyXvZTPvg9nZWVZXV4uaCgqr9U5k7OZa3Inutr2qrCiHStQLpa4rt5/rdvaFfUi6lmXl87RbXRSqqm7adVVuHzekW2gc3t7ezokTJ1ytX63T2FZFiHQ6zdTUFPF4nDNnzmzZrlutpvfkkUFujt0l6MuJ0E3T4GRPf9H2kiTRe+AwFw6cdbX+fkKtUgKyLOP1emloaODYsWP5nzvNHo4j2Pj4uCuntr0c6db6xlkLA/Pu7m6mpqby/5+enubQoUM7vi/sQ9J15t27QTU53c0aJIQQrKysMDo6itfr5eTJk8Tj8YpSGNuxdyyFwjTCgQMHCAQCrtp1q410P/m+n+X25D3uTo0ghOANpy/ykXe8Hz2r53WYezVCrQV2WjJWqtkDyju1eb3evHrCqWvUgnxrTbq7MYWiEtK9ePEiQ0NDjI2NcfjwYf7kT/6EP/7jP97xfWEfkm4lqKbDrNQ+zpSG0dFRgsEgZ86cycuL0ul0XmDvBoqiVDT3bLN23dHRUeLxeD6N4Dymbmfdcts6JNoYaOAr/+fnGJ+fQlYU2gPN3Ll1h3Q6nV9PCJHv9NpOrnIvkvcrVfzayqltdXWVdDrN9evX8+d/ffGukuOuZeS8Ww5jsVjMtZeuqqp8/vOf59FHH8WyLD72sY9x5swZ/tt/+28AfOITn2B+fp4LFy4Qi8WQZZnPfe5z3L59m8bGxpL7ukWddDfZp3BKQ2NjI/fff/+G3u7dNDKHjWTrSNGgcu+FagppkCuWHWo5yPDwMEvpOQYGBmhoaMhrpEdHR/PmQTMzM0XGNIUuYfvRyHwvNUcUOrVpmoYQghMnTmxwanMcwSpxaqsl6e6W5rfSjrR3v/vdvPvd7y762Sc+8Yn8vzs7O5menna9r1vsu6u+ooR1laRrGAazs7OMj4/T0tLCuXPnyhqH1zpdsB5ORLoZ2a7f1g0qae0tjDgL0xn9/f20t+f0uE707uQqVVUtqiSbpkkikSipbXXIOBQKFU0v2IvYqybmhSRZ2OxRmGqqxKkNatcRuFukaxhG0dSMvYp9R7qVoNKLxrZt4vE4q6urdHV18eCDD275IW7VkbYelZKuYySt6zr9/f0lydZBpaTrFrIsY5omt2/fJhqN0t/fv6kErVRaQFXVktpW5/HYcQlz5FTBYBBJkrAsi2w2+4p455bDXtTWulmrEqe2ZDLJ888/v0FJUQ151r10i7HvSHcnvni2bTM9Pc3U1BR+v5/u7m4GBgZc7btT6QUnsl1aWuLw4cMMDg7WtF3XLXRdZ3h4mFgsRm9vL6dOnaqJibmzXTkj82QyycLCAslkkjt37uSnSBRGxdWSwHawl9ILhdgOgZcq3l25coXTp08XObUlk8n8CPhKnNrKzTOrFpvliPfKjXkz7DvShcoLLOUu7sIpDZ2dnTz00EOEw2Gi0ajrtXdCjTA6OkosFqO/vx+fz0cgEHB1MVViYr4VHO+KpaUlent7iUajdHZ2utq3Uley9ZBlmYaGhnzKwrkBOnaNiUSC2dlZEonEho4vJ0WxU1++vUKUO7kWlHZqg+IR8G6d2nYjvZDJZPZFagH2KelWglI6V9M0mZycZHZ2lsOHD3Pp0qX87yvNA1e6fTnSXU+2ThohlUpVVBzbbouxaZqMj4+zsLCQt54EinSJW6GWJuaF0DSNlpaWDRX8wo6vhYWFfNGo8LG4VjejV2Okux6bHVc1Tm1OS/7q6mq+2WM779s0zQ0F7Wg0ui8cxuA1QLoOKToFsomJCRYWFjh8+HBJ4/BKSbRSg+9yaoT1ZFu4/U5M+F0P50Y0NzeXN1V3vsSOKqES7FaOrVzHl+MQ5lgCZjIZrl69iqZpRSmKrUxp1mOvEuUrrTbYzKltdHQU0zSLrBmdVFGhT67bgZKl0hX7xUsX9inpVhJJqapKOp1mcnKSpaWlDYRSavudMqSBl0l3K7It3N5t80U1ZOCkWKanp+nu7uby5csbLuhK163lo3215L3eISwcDvPQQw/lTWkSiUSRKU2hrnUrE/O9GunWqlhVSwJXFAVN0/Jz9hwYhlH0dJJMJjFN05VTW6mJv7FYrB7p7gU4lfEXX3yRvr6+TcnWQaWFsUqRyWRIp9M8//zzW6oRYGeKY5D7Yum6ztNPP01nZyeXL1+u2Zd2uzndwnVqBWetcnnKwkm3zqOx44NQqGvdD5Kx7aLWOdhS62maVrVTm5O2KEQkEqmT7k5iqws1nU4zNjZGJBLB5/Nx5MgR15Nsq9H2ukFhZKuqKpcvX3ZdHKtl5O0YrUxMTGDbNhcvXqzphADYHxXkQkiSVFLXWkpKlU6n812FhSmKagjv1ZRe2Go9Nzd0t05t8XicW7du5Zs9rl27xvj4eH5bN9ffVibkQgh++Zd/mSeeeIJAIMAf/uEfcv78eQCOHj1KQ0MDiqKgqirXrl2r4GzsU9ItB2dseDwe59ixY5w6dYqxsbGqmhEqwWb2i6WaGp566inXxFSL4hjkLqK5uTnGx8dpa2vj4sWLXL9+vabRbeG/azl91w2SeprR8DSKJDPY1oumbP99lZJSzczMoOs6jY2N+Xzx+uq9Q8abtT+/tDjGr/797xK1Upw9dJJPv+1jtAWqj9T2KoHD9iVj6/P2q6urnD17Nl9ovnXrFi+99BITExN8//vfJxgM8o1vfKOo0FcINybk3/zmNxkaGmJoaIgrV67wC7/wC0X2jT/84Q/Lrr8V9iXprr+QE4kEo6OjpNNpjh07ViTcrzRyrSZKc6LRwgvVTQeZ27W386jutDKPjo7S0tJS1PBRiW1kJdgp9UI5LCfD/Mb//j1WU1EEgv7WI3zmLR/Hp+6MhMjj8dDW1lYUjdm2ne/2Kmx/VlW1iIiDwSBRPclnvvU5ktkkLcFmnp+7y7/77hf5nfd+Zk/46e5GemE7cN6rJEmEQiE+9KEPMT09zU/91E/xkY98hHg8vqGgVwg3JuSPP/44P/MzP4MkSVy+fJlIJMLc3Jxrv97NsC9J10E8Hs/bK/b19ZWchVaNkXmlcIpvmqa5Jlu3j0HVphcck56RkRGampo4f/78hlbmncoX77ZRzZ+8+E1W01EOhFoRQnBvZYLvj1zhPSfeUvPXKjcBwTH3cVpoHRQWjObn50kkEtxcGSESj9Kg+rEti1ZfA0PL4yT1NCFvYMPabmDbds2Iba+TLpQ2MHfyw+s77tbDjQl5qW1mZmbo6upCkiQeeeQRJEni53/+5/ln/+yfVXTs+5J0k8kkt27dwrZt+vr6yk5pgFzCPpFIVPwalRQ5FEUhkUgwMjKypRrB2d5tnqtSYnQmEYyMjBAMBjl79ix+v78ma7t5bWd6RzWkm9RT3FkaBeDUgX7X6ywkVghoufcoSRKaorKYXKn49d2g0vdVqmDkX2jmq5PfAyuXloolE2SyGW4+/wKNDY1FUbFbTeteTy/stLlRLBarqYH5Ztv8/d//PYcOHWJxcZGHH36YkydP8pa3uL/B70vSVRSFvr4+Vye5msKYQ0Zu7s7pdJpYLEY8Huf48eOu0giVkG4lka5j7zczM8PrXve6TR+xoHKNsVsUkqUQgunYAllT51BjBwGttHHQajrKr/3NF1lNRQDQZI2Hey7TYgc4wUaD+MXECr935X8wG1tCliQieoyA5sMWNoZl0tPUxX9+8o94cX6IroZ23hxwN7/K7fvbDk519PH2/od44tbf4pO9yIrCpx/+OBcHL+ZTFKurq3lNa6GBebn2572cXqi1zK4UKmmOcGNCvtk2zt8dHR28//3v5+rVq69+0vX7/a7vnNuxd9zswitUIwSDQXp6elwn1ishUjfRaCQSYXh4OJ8/PH36tCuhea2kXU4Xm2EYNDQ0YNt27o+w+dL1P+fvJp9DliQavSH+zZs/TmfDxqnAj9/5AaupKB3BNkZXp5mIjnBvcRy/5CPhN3nw8Gkypk5vcxe2Lfj5v/h1FhNhvIpG2szS2dDGgrSCJEm898Tb+d7Q0zw7+xKN3iAvzN/jRvYOb7r0Rhq9m9+IACLpGH/wzJ8zsjrFYFsvP/fQT+b3KySQpJHm3uo4trAZaOmhybv5Y60DSZL4v972MXpoI9jeyMmDfZzsyE2PKGVIU2hgXuiB4Pf78yRcy6nXu5EOqBab2Tq69dJ1Y0L+3ve+l89//vN8+MMf5sqVKzQ1NdHV1ZU/9w0NDSSTSb7zne/wq7/6qxW9h31JupXaO1Y7sqdUL3eppobh4eGKyKsS0t1s21gsxtDQEJIkcfz4cRobG7l+/XpJ8XgpVJO6KDz3hd4Vhw4dwufzkUwmCYfDJJNJrs3e4lszT9IRaMOjakQzCb707F/wb9/68xvWXk6F8Soe0kaGmfgiHlnDo2qEZD//6ck/ZOBADx5VpcEb5B29l1lJRfPV/oDHx2JilW/+49/Fp3rRTYM/vPY4BwLN2JJNSA2wkF7i3tI4F7o3N5s2LJN/8+3/ytjqNAGPn/HwLOPhWT73E59GkZX8OYhlE3zhua8RzkQBiaDm55fO/RTtAXdffFmSOdNyjNMnTm8p2StnYF4oo0okErzwwgtF2uJKZqwVopaNFrVGLUjXjYH5u9/9bp544gkGBgjCFKcAACAASURBVAYIBAJ8+ctfBmBhYYH3v//9+WP56Z/+ad71rndV9B725pmtIWo1smezDrKdtHcsRYyJRIKhoSEsy2JgYKAozVLt7LOt4KQMnOjY0fp2dXVx6dIlIHcRdnR00NbWxuzsLDFsPIseZFkik80gsia3pu7xwgsvbPDQvb/zOM/N38EWXhACW4KgFiCjZ0mZaRq9QYJeP6vpKN8d+VHRsTlPnDcXh8mYWQ41HECSJJIiRVrKgAARhLnU0pbvcyo6z2RkjrZAc85nQPUysjrJbHyJI02d+XPwd9PPEsnE6Azmnm6WUxG+N/EUHz7l3th6O4/d62VU4XCY+++/P/e+C8xonCcQj8dT1OQRDAbLpiP28iTgck+g2Wy2bO2iFLYyMJckiS984Qsb9uvr6+PGjRsVHPFGvOpJt5rqf2FXmpt23Z2cHlH4WslkkuHhYbLZLIODgyXv7DtFus62i4uLjI6O0t7ezkMPPZSPqAvfjyRJWLaFoiqYWKgeDZ/PR0Y2efDgKQYHB4sMzdPpNK2yl4tNJ3lq6UWEEIS0AA1aiPHENJqs4ffknjqCmh/LFnQ1HGA6Oo8qqxiWyaGWA/w/z/4ZAglZkjjbfYKnF28grR1Wg+rn2tJN3tX/pk01vJqsYhfkDQU5UtdklaSeJmWkCUpBEkaqaB2vqhHTK5s8vRPFL1mWi9qf4eVOL4eMnU4vIUQ+ReEQsc/n2yB/3O5x1bJZZr976cI+Jd1KPsRqPnDHr+HWrVuu1AiqqlY096waI/MXX3yRVCrFwMAAra2tNTMyd3PBCiEwTZNnnnmGpqamLc3ds5bOf3ruq4wn58naOjPRRfrauultOsTHzn8g71JVaFBjmibHk8f5yfgjvDQ/ypduPs5cZIEDWhOSVyaTyeDRPESzCR46fB8feedP8P8++1fMxhbpbGznxvJLdARz5yVr6synlzl14Bi2ZRPQ/Hh0GVNYJM00zUpxzjRj6nx/7Clm44sca+7mgc7jfHfkaQzbRAibgbYefuOHv8ethREMw+BtRy/wwQcf4Zm5m2RNHUmSiOspHmkfdHXeC8/rbvg4FHZ6FSp9Cp3BYrFYfgy8ruukUinS6XSejN2a0azHbna37ZdOyH1JurBzWtB0Os3y8jKGYXDy5EnXaoSMkWVodQIbm6NNh/Eq5fNobkk3k8kwMjJCOp3mxIkT+dE4Wx1LLSPdSCTC0NAQ2WyWBx980FWF+C/ufJ+XVsc52JQj1dV0lFOtA3zmrR9DkUt/AQu7v450H+HhC28lGo0yMTHB85lRnhh5Esu26PK0cr/oYWFylg8cfSuhUIiXYhPcXL33sr+ComELQXuwhYDqw6/6GJuboNPTQkgr1sFats1/vfoV7iyNspqKspqKYQsbv9eDntGRZJiKzXFnaZTDDQdp8AT535PXOXVogPcffyffn3ga2xb8eN+buXTofjenPY9aT8mtdK1yzmC3bt2itbUV27Y3mNEUpijctD/vxlBKZ0T9fsG+Jd1KsdUFXphGaGzMaSXd+jVkbZ0/GvlLUhM6CDgYbOdfXPg/aPCUrpRvRbrZbJbR0VHC4TD9/f1EIpGiqHAzVDqcshzpxuNxhoaGADh58iT37t1zbRI9Gp5BkZT8+fapXuZii3nCtWybucQitrA51NCBKqtrP7dImRn8qg9VVpBlGVVV+ZmH3sc/OvfjGLZBgydYZNs4Pz9PKhwhlUgxk84S9ARIWike6DzJe0+9na+99ARxPYkmqXz0zD9EXUf6U7F5hlYmCCdjzERzOd+0kcGnegn5/Pg0D/PRVRRJIZyJ0uVvxyNr3Ji7y0fOvoc3HD7n6pyUw16MzoQQNDY2FpGxM1qp0C6zsP25MEdfqC3eDQPzaDRalE7Z69i3pFtJpLuZLrZUztaZnuoWP1q8wWJmlWMdvQDMJ5f5zujf88GTj2x6POuh6zpjY2OsrKxw9OhRTp48iSRJjIyMuD6W7eZ0U6lUUd7YKdJVcr6Pt/Xyv0eeyW+vWzrHDxwFIGvq/M61r/LS8ii6ZWAjOHfwJKcP9HMrOkzSTOORNT7Y/wi3l4a5Mv88HfozPNLzRs605iZIrPdEOMEJjhw/yn9/9i9YToY53dTHm4P3ER5e5F2eC6h+D0kpTsD0kNEzXJl7kcXkKn0t3TR7c1/WmdgSXlXDEjayJWPaJootIwSosoxprr13AYZt0d100NW52I8opVEvHK1Urv25cNil45frBAG1injLGZjvFy9d2MekWwkKjcwdbKVGqKQwtpqNohacSp/qZSm1WnZ7x6XKgWEYjI+Ps7i4yNGjRxkcHCx6bNvMUGc9qiXdTCaTPx8DAwMbWqoraaT44Okf46nh5xhJzQAS93UO8k8vfACA7409xa2lYRo8QaZi8xiWQdrM8KPl5znaeJi+lm7SZobPv/hV/LKXgOzBEhZ/Nvwtmk7/JN2h0iODBlp7+K0f+5c8P3+HH05c4YrxEu8aeDOd/jaemXyR0cQUgakQXx/5DsPJaVRZQZJlHum5TLu/BdO2UCQZGxtNyalRBIK0kaGrqZ2FaE7SFjeSdDd18tHz73V1LvYjKolOt2p/dgqlN27cwLKs/Hw1J03hdhSVg3KR7n6xdYR9TLrVaHV9Pp8rNUKlErCBll6uTN3AFjkCSxlpBluPlt3eiXRN02RiYoL5+flNzdUrMaaZySwyHH6W9ng7bzp8njZ/+QhAkiRM0+Tu3busrKzQ19dXdvBkJY0UXs3DP3/dT9J+rAtrLSqUpdyxT8cX8MgakUyMrG0gS2BYBj7Jy3I6F336VR/RbILWUCOyJeFXfSSMFOOxGTySxkvL4/g1L2cPnkBTXs7lPTd/h9++9v/hWVMgXJ+7hWnYxLNJ0uk0341dR9VkjrZ3Y9s2WSPL96au8ImjH2AsNMlQfAa/6qHJE0JVVS51v46onqA92MpDh15Hq6+Z6elpLvefpcnnrhFiP6IW6gWn/VnXdTweD8eOHcvPV3P8chcXF0mn03n52/oURSmUIt1IJFKPdPcaNE0jmUwyOTnp2huhkkj3zUcu8MLoLSZTywC88fA53tpzsez2kiSxurrKwsJC2WkN64/HTdvwreVhvjb5LRRZQYtrXJm9wb+++LMlidc0TVZXV/OR7froutQxu410nfN6uHFjTry38RB/N32dqcQ8GTNnRJSxdDq8rXgVDxkzy3h0hng2waoapV3KRTC2EMzHV/jiM3+CaZtgSxxvPcqvvuUX8KwR73dG/x6vrNHgDZLQU9xZHMUwTXqbDiGrEM7EEBmbztABFEXBL/tJC51z9z/AHz94mW/c/gHPTN2kWQ3yjoMPohg5AvL7/YSUECFviECjhE/bWwMQCz8X0zaL8unVrlfLlmLnui2cr7Z+tJKToigc6eOMVirUFpeLdN02RuwFvOpJN51OE4lEWFlZ4cSJE67UCJWmFzyaxpuaHuDcxZzJcTnlgm3bTE1NMTExgaZpJWe0lYJbtcO3x57Ep3gJqD6CgSCLqRWuL9zikaNv3HAM09PTBINBent76e7u3nLtrdIWbvx0LdvC8hhIfptmJUBW14jFc3rRcDxGT1MXV+dexLRMmn1NTCUWiKtJMnGTpWSYZ2ZfRJJBU1VMy+LJ6et87eZf89EH3pc7RklCIIhnU9xeHCFj5FI40/F5WuQGgj4fcT3FaiZKQPWRlOK0tDTwvdkneXPXJT74uof54OseLjpmp/PLMc+ORCJEo9F8s0Hhn1eqoUAIQdJK85s/+j2GVicIan4+/sBPcq7zVNXr7WZLsaIoZduf14+ATyaTWJaVL3Y73+/9VEirrZHqLsLN9Ihbt27x/PPPEwqF6O/vp6Ojw9XFVM1wStu28SqekoTrEN1TTz2FaZrcd999JU1Ltlp/K5jCQpYKc68Spp0jayEEMzMzPPXUUxiGwaVLl1yrM6DySLfUti+G73EvOkZQ8yNsCPr8NAeDqLKKJqt8oP8RWr1NvK79OP1NRxhs7EWyJZKpNKvpCEIW6KZJUk/ncq9ema/f/hZj4WkAHu17E0vJVV5cuEfGzCLLEjY2WcsgYiYwbItHB9+AgcFwfIJ4NsXRhm5urN7hP7/wB/z1xA/RrZdbxoUQDEcnuZMYwwyQd7Q7c+YMDzzwAAcPHkQIwej0GF948it8+juP8e9++Dn+8sZ3mV/MFWPdnjMhBEupcC4ar1AKads2f7X0d4yEJ2nz5zrpvvDs15hPbN2Bt9PYTgHNGa105MgRTp06xYULFwgEAvT39xMMBonH4/zGb/wGn/3sZ/n85z/Pz/3cz/E7v/M7m7b9f+tb3+LEiRMMDAzw2GOPbfi9EIJPfvKTDAwMcP/99/Pss8+63tctXnWRbmHO1vGznZqa2tGx6uWIXAjB7Ows4+PjdHR05Du4KhmrDu4j3TcdfpCvrHwjly/O2KiSwv0HjheZmBeO59npluH1WEyv4FU8HAy0MxtfQtg2Ho8HOZWlv7mHgObDq3jy6QKvooEEU9l5TNVCFWsFLpGTnSmKgmorXJu7SU/zIZ4Y/Vsyto4tBLYQyLKELMlYtkVSpGn3NvP92adJZNJIEqSMLH9190lOd/eiKQpPL94gY+l8sO9RhBB8/d43eXL2GrIk41U0Pnz8PQRF7tH95sowX735V6SMNKpfJSuyhO0YVkYwOj/LSHKKtwfPo2d1FEXZNCpOGmn+85UvMxyeRAjBm7rP8/Gz/wjF5SO+YRpMZ5bobs15vfpUL2kjy3h0ls6QO6nhTmEnWoqdc9jR0cHv//7v82u/9ms89NBDHDt2jJs3b5Yl+e1MjHCzr1vsW9JdT3SlyNbZxjEXd4vtWh4KIZifn2dsbCw/GqewMFBpR5pbwnv9obPEIlGuLrzIweYDPNTyOmZuTxAMBjl37lxJE/PtRq/ltrVsi6n4fE6LG+pAk1UO+Fq4HRlioKWHlVSEmdQC6WwGWVbwa15M06TRE2I1HUFTNOLZBIZsYlgGApBkCc2joWd1ZFnGY3vQFBWPovGlF/8HP5y9iqIpBIM+0noWRZXBhkxWx+v1kJazpFIZJEnKpyJ00yCeStPe2Eirr5lnl29zyH+QmJ7g60PfxBa5NtYmTwN/OvQt/knLo4zHZvhP17+MR1YQEqzEImiKSoMnhCxJZMwss2KF1r6D9DR0YZpmvu15bm6ORCKBbdtkMhnGxsb4X/N/x93lMVr9TQgEfzt1ncHWo7zj6CVX5xsBfsVLxtLxq16EEAhsGr2hrfct8/nVCrvh4xCLxejs7OTixYtcvFi+lrKdiRHj4+Nb7usW+5Z0HWxGtg52atjkeggh8t4Ezc3NZdtlKyVdt9tLksSDHWc4bLeRzWaRwuamvrqVKBIq2dawTf5y9Umi19NISHQFD/DPH/gI97WeYDo1z2hsisH2HtpSzdyZGaexJUQsm+B3n/0TfuWhn2EoOklMTxCQvfz55PdyXXZryhBFUVAVlWxWJ5nN4FU1fjD3NPOpFRCgyQpCUxGywLJsZE0i5AmALQBB7tIQQO4aEQiSZpp20cRCYpm74XFuLg5h2Ca2LVAkBRCsZCKokkK20WBoZQLLtgh4Q+jCQBEyumUgra0nSRIyMoade8xVVXWDkbllWVy7do1gMMjQ6gSYgng8nlOU2AY3pm5zofWUqzSUEIJ/2PVWngg/RcbMIoTNpa77OdXW5+rzKkStJ4nsBum61eluZ2KEm33dYt+Srq7r3L59e1OydbDTpOt4E1y5coWGhoZNpzXAzkW6zqy4TCbDuXPntiwuVJpe2CzSzWQy+YLG38xeZVpf4khDzux5JrHIX439kJ868R7ec+TtRPQYQsCvP/lFGj2hXB5cgZV0hNHIDP/45I8D8IPxp7ClXJrElmRM20JCEPR4afc0EfMkkVUZXZhkzCyqopC1dSxpTV63FtFato1pWdhCoGkqum5iIxAil6KwJcFUbIlENoW1Ru65aNHJk0sI2wIJfLIHr+TBlmxMTCQhIYtcGiNpppEliZAWpNEb4lCwI79WxtLxyGpRG7SqqnR0dHCys5+nZp6n0deYi4DTBgd9rUXeueslVYVDL23b5kRTL6+//yIT0VkaPUFOtvVVFbHWcuwPbH8oZSHKmee4VS9sZ2KEm33dYt+SrqIoHDhwwLUaoVJPXXBXxV1ZWWF4eBjTNDl79qyrKmqlH9ZWJJ1OpxkeHiaVStHV1UUikXB1HLXI0xqGwejoKCsrKzQ2NjI1NcX15ReQbEinUiiqiiapzCVzRR1Zkmn15qISv+ojlk1S+CzgV3NpmNVMlO9OP4VuGNiKjSIrqJICWHhkJT+E0jIEM8l5UEFRJTQ8pLNZPLJGZ/AA49EZJEBCQrZlbGHj83mwTAtJUfBpXvSsSWMoREysjXUSL9tFAnnXMUWWmc4uMi7PYHtNZoxFbCGQBBxr7KbBE8CwLAZbjvKBvh/L6Y31OF8d+itmUouoksL7jv4Y59pPFRHIT515D2PRaRaTqwgEZztP8Y/OvyfvYlaooIhGo3ljGqfry5E4HvC15K0mq8UrNX7d7Vqlji0Wi7ki3e1MjNB1fct93WLfkq4TJbhBNZ66W2ljw+Eww8PDaJrGmTNnuHv3bsVm0ZUey3o4Hg2RSISBgQHa29uJx+NEo1FX626nkGZZFpOTk7wwfRtPm4+e40foCx1BkRUiUzpDN6ZQ1p4wotkYPXY7165dIxQK0dDQQCgU4h+f/HH+yzN/xEomN6KnLdDMm49cAOAPbv4pMSOJD0+uXdiy0TQVTZYJJ5JIkiMRg7SuE/B7sYWNkHLRrCIUNEVBk1VUWclpgmUJW8i8ufM892LjhLNRJEC3DObjy0hyLs+bi3HXvX8pl6b40/D30FQVIUtIioyCwCd50G2d8+0X+cnBR4v2+9ORbzKbWqLV00jW0vniC39Md6CLQ4EDnLQPA9Dia+Q33/ovmYrNo8gyPY1d+WYSKPbOLbzmDcMgkUiwtLREOp3mueee2zIq3gq1tHV01qsViZdTQqRSKVdeutuZGHHgwIEt93WLfUu6lXakVWtkXkqIPTQ0hKIonDx5Mq8trHZqrxusJzynbXhpaYljx47lPRpKbbvVupUW0goVGasNSW74hpGScHX0JieajvG+o4/wlsMXePrudWZFGEmB84fP8E/PfBAVJW+YMjs7SzoR4YPNb2HaWqLBF+JNPefxopExs0wnFmjxNCB7BTEpRSQbQwIUZJDWjodc9GmaJpIkoSpqbhtZRhYyJiZNnkbi2QSWZOeOXxK8EL1N1silCyxhIym5FETI6ye1bnK0LEsoyMiSjCmZ6CJDTDcxbRvWyDmLQVRPcG3pZp50hRAMRUd5KTJMkydX0JqMzrOQWMEybWbi8zxv3uZB41zOelLR6G85QiXQNI2WlpZ8M8PAwABCiLKj4AuJuFyuuNbphVp6Bm/mpevmNbYzMaLcvtVg35IuuK+ob9fIHF523RJCMDg4uKHXeyfzxoqikM1m89Hl7OwsPT09XL58ecPFVmnKoJLCiWO12NLSwrkHz/O5W1+i2duIKqsI2+ZuZJS51BKHgwd5R8ODnDp7BgubZk9D/qbQ1NSE5RP8IP0UkWAMLaTSIbfywuo9bo8Nc2yii/vkQYyMTlSPIdsSIc1PTE0gJMjaJpYtQM6Rq7EWmSnyy1OIZVlCAsJWlLStE88mkWQJr9eDosjo2Jgilx/2ej1YQiALGZDobGwjnIrlfHh9HhRFAgGKpJI2U6AIbCNXlHPKcUKAkARpK5U/Vz9avMa15RsIyWRZXyFixJhLLKHICiFPAK/kYSm7wlBkkgcObBy8WQkKiU2SpLxd48GDL5vyOFFx4Zy19SbmoVCopjlYBzs9H62S16h2YkS5favBviZdt6jWyNyyLBKJBMPDwxiGwcDAQNncUTXE7rbzR5IkwuEws7OzHDp0aNO24Z3Q3kYiEaampvB4PPkiYcbMYgl7rbrPmgxLRrf1/P9LSZaEEPz11PdIGElaPE3MJRe5Hn2J/sZePLLGTHaZwUMD/Fzvh/jSrT8naaSJ2HEkJffojwR+n0Y6a2Cvne+A35NPNchOcYlcVK5qMqqm4vd5EIBHVRFCoPk8dDSEkOXc9qupNBnTQhIS/a3dxIwYPsVD0kjnZGOyhGmbSIgcwdsStnCIN1eU623KmfHols6zyy/SpDXSExK8FBkhIzIYtoFP8aLJ2prDmp0/3u3ATTTpRMXr56ytj4pTqVQ+37xVVLzbKEW6O3GT2Gm8Jki3GgghuHfvHrZt5123NkM1c9K2epRz9L6Ol23heJyt1nWDrUg3mUzmZ7F1dXXh9XpRvSrfGPsON5ZfYim9QkJPcjjUSdpI41U8dPhyhZxyN5OMlSWix2jx5J4UUmYWCbCFhSx5CagBhqLjfOL0TxNQvXz97uPo2SwZy8jZBNoCRZVoUjw0KQHSIgsyWFhYtsAUNrYAOaffQlMUTnR2MBdPoKm5c2MBBxsbEAhM20aRJA4EA8yv5Yq7G1oZisbR7QxIAlVWkYWEgYmNhKrI2FIu9WDaAo8kEfJ6ONuW02watolu6yxlFlg1wgQ9KoZt0+jzEc1kuBMeQZEkvLKH2cwC93GcmJ7gD+/8T4ajk7T7WvjoqffR2+CuUFPtI3ypqHh5eZlwOEx7e/uWUXElueJaoBTpOv7X+wn7mnQrnR7hJrJMp9OMjIywsrLC4cOHGRgYcHVhVTsnrRTpCiFYXl5meHiY5uZmTpw4QTgcdj3htxIT81LnL5vNMjw8TDweZ3BwkLa2NmZmZjAMg29PPsnVhRs0exo5GGhnJrlAMBvjcPAgDx96M0Ft84KGR9HQZJWspeNVPLnUhBCoUu5S1C2dJk8uTz6UHKbJF6CJIJl0NJc6WMvntnpDnGntZyUdJ64nmUktENA0vKqKJWyimSwSggavF1uGxoAX3bKQkFGRUNcIU5FyaQkkCZ+qYpJlJbOIYZkIkfu9YZu5gptt5qJpJCRZRgFCmoyQJGxJ8Gdj/4sbK7fRZIXl7BKWMDFF7rOQAFVV8Gkqli3QFI0e3wF+tPAsxxuP8fV7TzAWm6FBC7CQWuazz/0Rv375X9Do2brBoZZ5U2eS9PqouHC0j9tcca01v+VqLPvJYQz2OelWgq3UCI6fbDQazfd2Fzrgu12/kuMxTXOD4iEcDjM0NITf788/ykcikR0dNunANE3GxsZYWlraoH12bnC3I0M0aCEUWSEg+2nzNfOWQ5d466FLrubEKZLCo4ffxhPTPyBuJLGETrO3gYgeR5UVgqqfR7rfDEDEiOFBo83fhGlbRI0UIBFQPLyx80H+Qe87WclE+OyLXybg0fAqCpYQqLJMW8DPSiqNpiggwKdqeGQZUxFosoIqyzlCXSNWz9o+thCsZJPYwsYGECIXiRcpGnI/c8haliQkZHRhcTtyD1VRafB48vleCcisFe58Xg2v7M3pesnQQCNzqQXG47M0eUK5x3o5QMJIMRmf43VtW89dq/WAy1KBQOFoH7e54kAgkJ84Ucl3qRxM09zQVbnfvHRhn5NuNZ6660lX13VGR0dZXV0t8pOdnp6u2K+hEi3w+jRALBZjaGgIWZY5ffp0kSl0tdODLWGxlFnCsA2aPE00ao0btrVtu8h57MiRI5sW6EJagKV0OO+PIIQguG7umPPzcjjW0MP7et7Fn47/JV7VgyTJTCYW8MleFEkmZsRp8zUTULzMKBEaNQ1DBPAbGi3eZvoaevkHR36M33v6z7i3PMGitcLhw6E1RUFOV6vIMj5VBQlUWUGRJWIZM1eMExaqKdPg9eBfywXrprkWSUtokoSiqmTMnJG5T1PxyDLLmQymbecLaM4fpwEjR865HwoBkkxOWywBGNjkdMqOHMw5R23enALBFDaapCCEyOmJVXcWkrZt77gWthw2i4rD4TCWZfHSSy/lo2JHLujMWKvktUrlb+uku4exXqtrGEY+qjt27BgnTpwoIqxqJvxmMhnX2ztpACdvappmSVWEs22lj2qWsLi2fI3V7AoAEjLn287T4X9Z5ylJEtlslqeeeoqDBw9y6dKlTavDQgh+ovcdfOmlP2M1GwUEXcEO7m87WdGxAYzGJxBC0O5t47nlu2v5Vw0B/NG9v+BnT3yAtJ1EJdeG2+T1crKpj8sHLtGghfjn3/t3JIwUS/Eo0WSaQ10n1tpwc40QEhByjH0kSOo6NqCsvQ/dtlDll8lPkiScT1+WZAQ2qiwhywoBTSOayeQiSnLEKoTAu9b15kTDErnPyqvIa8chg2TT7GlEw4NtJ/ApXsJ6HEtYNCohTjUf428Xn6Q9FGQkPEdQDaBKKucOnKK/0Z2ErJZa2Fo0MzhRMeRI0ZFWlRsDv9mMtUKYprkhxVZPL+wyqtHqmqbJ+Pg4CwsL9Pb2lp3WUG2O1i2EEPlOtq0KddUoI5Yzy6xmV2jUGpEkCd3SuR29nSfdlZWVfJHs8uXLWw6ddEi3p+Ew/+K+jzIem8ajqJxs7sejeDZEtlt9NlkriyLl/AkM20BTNGxh4Vd9xI0kw7ERNFnloHSA5qYWTMsgoAXoDnbx+Vv/HTQbr6RwtKuDpXAMXTdpaQhgC4Gx1vLr92g5QgWstc4xR+PrkRWEAEXORboCUMgRqqN8sIGQquZuTpaVz/8KXnY5C3g8+UhXt3Itw7YtsGQb3TDxKz50yyBlp/F7JbJmCq+ioEleVCSmkhO5G9eBAVq8ITKmxY93v4OzHaUneJTCbqQXqsH6HKzH48Hj8ZSMiuPxeNGMNcfAvDBXXOrmUifdPQxZlpmeniYajdLT01OWbB1Uqrt1u72TzlheXqa7u9tVoa66SNfMRXxra6uyStpKE4vFuHfvHqqqcv/99/P888+7mvJbeAztvhbafdtz6h9oPMYL4Ts5GZYEmizR7AlhCRtbWIS0IPPp3Otpsoph63hlL1PJGXRLRzdMLNvGsnVOHD5IUzCQy98C3dFkMgAAIABJREFUKcMgls1imCY+RcWrqaRNk6zIfT6Or4IiSQhbYEs5wrXWolZZkkCW8ctyPgKWpJxErOiKkSSCBZGXLQQpwwAhoVu510rbGdJ6OucfIXLkriq5YpyFyVxqiZ6GXGfakcZOonqUM+0DKJJ7Eq11Ia2WUfNWa5UbA+8YmBdGxU6+2DEwt22bcDjsyoR/L2Hfmpi7hWVZjI+PMzs7i23bvP71r6enp2fLi7RS0t0qGjVNk+HhYZ555hkaGho4cuQIDQ0NrqKZaiLdJk8zsiSTsTJYtkU4HcZasbh79y6Dg4Ocuf8Ms+Yc8/ICy5nlLderRCmSMBIk5SRpq7yd5tGGI7yn+50EPX7OtB/maHMrrUEPspzh7V2XeaD1PkJqiBRpInoEW9jc33ofsiQT8PhzvgQSeFWNBr8vf2y2EHgUJR/dJk0Dyxa0+XOkLNaIT5JykywylknGNLEgT7iKLONRVDxqrsPNtCya125Mhm3nZWZNXm+ewB0S96kqQsoV57yKioyc93Gw1v7xcgpEImPprKRXGYqNcS86QkxP8MzSNa6vPEfadGdHWuvxOrtJuuXgGJj39PRw+vTpvIH50aNH8Xq9hMNhPvnJT/L7v//7fPazn+Vf/at/xVe+8hVX1+jq6ioPP/wwg4ODPPzww4TD4ZLblTMt//f//t9z+PBhzp49y9mzZ3niiScqem/7OtLdjLBs22Z6eprJyUm6urro6+vLdS65vAiqiXRLEaNt20xOTjIzM0N3d3c+wh4fH69I2lVppBtUg1xov8jNlRdZCC/iSWqcO3qegwcOYtom357+DsvZJaJKjCemvslbOt/C0YbebR/DcHyYO9HbhD0RkvMJLrRe5KC/9LjyE80DzGQmmM0YYIMhLFq8Cq9r68cje3iw9RzXVm1OtpzgUPAQDWoDd6MjhPUwh9paSGbTmGvn0LuWBnAKWposY4lc0StjGjT7/BwKhTBtx81GkLKsYpeFtfSBsnZdWWsphlwBTKLV7ydtGAgh8GkaqqzkrByBgObFI+e+ThlDJ20aGLa1Zkkp1qRjL7+aLnQUFHyKl8nULIqkYAuLiB7nxupNQh4/Y4kx/kH3T5Qd/+TA8UsQQhDWI1jCptXTXORo5ha1btutdeNC4Vifr33ta3z605/mkUceobGxkdHRUVdBzGOPPcY73/lOPvOZz/DYY4/x2GOP8R//438s2mYr0/Jf+ZVf4VOf+lRV72Ffk24p2LbN7OwsExMTdHR0cOnSJTRNY25uriIj82qaHQpJuvA4Ojs7NxSpqlUkuIVpmoSnwyjzKm849ga6TnXl15lJzbCSXabZ04wudPyKj2vLz2xKum4i3YSR4Hb0NkE1QIoUmqRxPXydR32PEjfiDMeHsYVNT7CHTn+ue2sxu4AkJDIim2tWsHSmk1PcidxhLj1PUk6gJlR6Qr3cjQ7xo8UrNHkaWc4s4/dqSHjwa16QXo5Ss6aJLMvIyOiWRa6kJa9FoDaGEFhrPhKspQ8K35nzu8KinE2u0cKrqvm8sWFbeFUVj5z7syYyw6OpqIqCJWzCmRS2nVtDkSSstbUkwCd5CNo+FEPGr3mJiBgKKgkzQ1ewg8XMEk8vP82JxuP581UKtp0z+vlfU99mLDGBhESbr4X39/wEfnVrI5hC7LWi3FaIRqMcPXqUc+fO8c53vtPVPo8//jh/8zd/A8BHP/pR3va2t20gXTeG59ViX5NuIRkJIfIO76WmNVQj6aom0hVC5MfjtLW1le0ic/wUag3btjFNk6effroosi6EYZvAy+dOcVy4NsFW5jiWZbEQXkAIG5mcj4Eqa2TMDCuZFX60/KO1qFBiKjXFG9rfQFcg56YVs+KoUu5R3kAwmhzFtG1CSggTg4ge4RuT32AkPkXSSAECJIFX1rCFTdrU8akqkiyjyQqreipPbLkoNDf2R0MlYWcQ5Mb5OIS7HoJcmiJ/fUkSwrZfbq4RuWg5qet4ZRmvppCn7XyxLpe3DXm8uZvBWgHONCGhZzGESa+/k67GQ4zFJvFLXiLZWC7lkU4zYU2QktM8t3iDu+F7XGp/iAfa7y/7md/9/9l79+C47vO8/3Oue1/cAQIgQQIgQFCkSYk3y3KcpJ4qdVI1bTpJXLtt6jhKWk8zVetbpMpq5M7U6XTsTmecS6WOPRN32rTO/GInrS2riRPFF0kkRYkWRRLEhSCuxB17v5w9l98fi+/R2cUusAssFInWO5OJRezuuezZ5/t+3/d5nyc1wWTyDk1aMQtcza3xg6WXebj3b237vVb6HhsJurUM9OwlduMEvLS0RHd3NwDd3d0sLy9vec1OouW/+7u/y9e+9jXOnTvHl770pbrO4R0NukAJyLW0tFR1a9iN2WQ9026CfnXx4kWi0ShnzpzZQuQuP59MJlP17/WGcK2YnJzEcRzOnDlDMLiVPwvQ6e9EkRSyZhYLi2QhxXB0eNvPryaOIxa7qakpfBEfKTtNykyD6bAaXyGgBZiOT2M5FpFNQMhZOcZT43QHuzkY6CNmXMPZzFQD+DHsAqqjguRgYZMwNijYJlkzh429mXtKRYUwijzcoBpA31z4usISqXzebXT5FM1tbInGmCLL2JaFhIQoMhRVzIoaD/bm9QqfMgnJ1UmQZRmfJNOmBwhomuc5KU5C+BQNmWJ5Q4CvjYNtO8iyRcSnEVIjaJqJ6aRBggw5NEXHJs+BpnY2jHV86DSpUSzD4nuz3yd9O0U0GHU7+pFIBF3Xiw2lQgxV2qSw4aDLPlZya/U8QsDbtylXr4D53/7bf5vFxcUt//4f/sN/qOl424mWf+ITn+Cpp55CkiSeeuopPvWpT/HVr361ps+FdzjoWpZVM8jtRlO31ojFYoyPj1MoFNyC/05Rz7juTrGxscHY2BjhcJgzZ85w/fr1bX84UT3Cw70Pc3nlMjHiHG8a4Uz7mW2PUam8sL6+ztjYGE1NTZw7V9TB7cp1cWXjFVbX14jqAQY5yvTqNHEnTp4iFchRHfJOHtu2GYoMMZeZBtlGlmRUWcWybdaycZJ2AkM3MB0Tn+xHlTWwHDbbXshI+BWN5kDYFSr3az5URSHqK3bDbcdGkVQ0qUDBsciYeUAioGjosoNhmUiOg0JRDCfqD4EDmUKO9Vx6s/HmuPq9bB6nJRhGluQi6CvyJvfXJqT60RUdy7HxSxKqrJApZLHtok6vtQkeGTNHUA0Q0HwMNx+kVe/CcWx8isrd3AJZK0NXoAttc0xaMRVOHz4NBltGcXO5HLZskbVzaJKKLMtkrQyDkSN1P0uN1r/dby3dVCpVMkgk4i/+4i+qflZXVxd3796lu7ubu3fvVtTl3k7w3DuR9+u//us88sgjdV3LOxp0VVXlzJkzNYmH74f0opB7BBgZGeHatWs1AS40RpUslUoxNjaGJEmcOHHCffhqaXp1Bjr4u30/x0vzL3Gu/dyO2Y33M73HFXbyjuOQy+XoCnTxd/wf4sriFR449gCaptFr9PKXS39Z5L5aNkbBIJgMcuXKFSws6HCKHmOShO2AX9Hxq74iWDk2iiITVaOs5OL4FB8aFgXbJKBqBHUfqiy72a/j2KiyuqkA5uBTdGzHRld0dKlolZY2s9gUdRyafIEi+8Cx8Ws+17I+oOkcULTNDKt4DzKFPHnbRNmcKrNsC3mTESHLCmbBQPNpmHbxOTOdIlWvOMW2qXrmYVbkLYOAHCBtpfnpA+9zv9u8lef/m/kGBdtAlmSyZoZmvYWAGkTWpC2juJcvX2a4a5h4LM6d9ByO49CsROmx25mdnSESKWbHtdRXa1W+qyUaWdPdq5auN37+53+eP/zDP+Txxx/nD//wD/n7f//vb3nNdoLnArABvvGNb3Dy5Mm6jv+OBl0A3yZtZ6fYLehWeggzmQwTExPkcjmGhobqrilB/aArMk1Jksjlcq49T6Xj70Z/YacHV5IkLMvixo0bJJNJhoeH3eM6jlNs5jgOhc3uvmzLbtbUpDfx010/zURyAsuxOBw6THeg+NBej11HSckErAC2Y2MVLNJmGl/BR0gKkTGz5LQsBcegN9jBfGYZnxzgQDCCSR7HKZYI7E3ilk/RCch+2nztxRFoiqWKnGVgYRY1GFQVFQVFVtjIJ1DkYoYtSgkODqqkEFCVzYZa8fkKbILyei612Rgr1pDZLD0ENX1T6B1Mx0aRZFSpWB+2nGJDT+TKslQsj2TsDG2+Nhwcrq2PspRdpdPfxt/pfpgfrPyQuJGgO9DNB7rev60MZEd7Bz/X9jBZK4ft2ISUIDkjhxE3WVpaYnJyEsuySpTCIpHIFqWwRoPuW5E113u+jz/+OL/8y7/MV77yFfr6+vjjP/5jABYWFnj00Uf59re/va1o+Wc/+1muXr2KJEkcOXKEZ555pq7jv+NBt9bYrZC5d7XO5/NMTk6SSCTcKbJKxnb7wb0VY8azs7Osra1x9OhROjo6Kh6rkfKO8KY1TywW4+DBg64+hQBb8X5N00gkEkxMTBAIBNxrtCyLsBTmdOQ0iqK49WFZlokXNvDJARzFRnZkChQwbRNF8aFYCpKTRS/o+I0ACS1Br7+dVq0Vn+Yj42RImglMp4BD0bm3CJZBmrQmDLvASn4Z27YwMdEkDQsLySk6SOiKhl/TkABd1op1fBk3K/VGUXFMRlYgrOlYtr353qKOsOM4mxSt4ufKtoNPKYJ03rLQFMnlCIODjIwiyYTVMBfaLvAnd57jyurrxQEMSeJM23v4xf6/VzOgyJJcbNypQaTN8wn4A7QFwsjdb447C6WwZDLJwsLCFk0Ey7IaVtfdb6ueanXenaKtrY3vfve7W/69p6enhHNbTbT8v//3/173Mb3xjgfdWkn7uxUyNzeFUKamplhdXS0RxSmPnZTMKr22lrAsi1wux5UrV+jv72d4eHjb62mUZY/Xmqerq4toNMqBAwdwHMdlaohFRixIuVyO4eFhV+NUnIcAZwHU4roCchDbXkOTNUyKWfKR4AB3s4tYkoWlWhzvPM5ybgW7YKPZGhkzw3J+GUxQVAU0scV0CKh+FElmMbdIQPHTprcRN5LkTRNFUrAcC1kqGlRmreymyE5REMfYnFiTNsHRm+W6z9kml9e0DQom+DQd27GKmbYjBhUUCnYOGY1kIYUkbVLPpDe3xAoyFyIXeH/3Q8SMBK+uXiOgBIpjxI7Da2tv8MGe99NW4+SfJEn4nAA5srApJ+nDX+K1tm6sM5+fR/Np9Lf141eKPZBCoUAymSSVSlEoFLhy5QoAoVCoRKCmXibCdk4P9Ual31UymXQ5u++keMeD7n6GoihMTU2xvr7O4cOHK6pvlb++kaDrBT1Zljl9+nRND1k9TbpqrIS1tTXGxsZoaWnhwoULSJLE2tpaCeAKIPLKQTa3NZGxMqTMFCEl5N6v8vsmAHgkepy4GSdVSAIKXXonJ8OnGAndR7wQZ2J5gvuaTjCV/hZhLQybFjuqqaBKGjFzg2KOV+TlZgpZHBlkWyEvG5yJnGNNX+PS+kVkZHRZ32QsWEgS6JIPv+wjaSVRkItZomSDbVE2OuGqhwU2682mY6PbNoZtIrhiDg4tvhYsUyZppHEky5V4lIvoC8Bh5TAnIvcV3TasApIku+UDXVFA0slatQsoAaiyRshRsDepeV7AvZu9y18s/sXmsAZcj1/nkd5H8Ct+NE2jtbWV1tZWFhcXOX/+PLZtk06nSSaTrKysMDU15UoressTfr+/agKw3/5osVjsHae7APcA6Nabwday/RfTbGtrazva43ijnoGK7UDXK2Le0tLC+fPnuXnzZu1bzT1o6iaTScbGxlAUhdOnT7uNQdM0yeVyLC0tuVSlhYUFZmZm6O3t5cKFC6StFC+vvYhpmzg4HAr2MRQ5VvG8xY8xrIT56e6/RbKQwLEdwmoEHLANm9RSiqgSRbZlHMcm72RxNg0mJdnhvuhxRpOj5OwsmqRh2AaGncen+QjLEQzDYCI2TtNGMwesblaaljElE03W6A4cYMVYLTIQJGG6CT5VQULBkVVMyyRnGS5QCoqaIinImzZFBVuIlMtulpo1s/zswZ9hMbvM5bWLmyUGA3NTW7fN18T9yin3vrT7W2jSI2zkY/RG2ghoKpqislKYp9NuQ5dLG8UpM0miEEeVNNp9HSV/k6SiuHp5XFm/goJCSC2yOpJmksnkJCea3zRX9IKkLMsl01+A2ywVmgiLi4vkcjkURSnJiL1C5vvpj/ZOFLuBewB06wkx8FBtm+TlnHZ1ddHT00NHR0fNdal6Biqqbevj8ThjY2P4/X5XxFx89n6MDYvXiuZcOp3m2LFj7sPsbZINDw+zsbHB5OQk6XQan89HZ2cnfr+ffD7PG6nXcRwIqkU2w2xmhg5/Jy1667bnoEgKzfqbTbnxu7e4nZ7AH/TR1N6MrEoMhAa4nR0Hu5hTBtUghpVHk1QsScNyrM3SgYIsyfg1P6qi0qQ1c/LQe9Buapw/cJ6XNl6kYBZIplNojgoO5OQ8sqLgKjVuDk9IskRxire4UMtIbu5bHBf22qS/yQGOqM30hQ6TMTMoskKLP0zeLGA6FiHNT1+kFyWuuACnyiq/ceyjfHP226AYaJLO0egR0laKVzYu0aZ10BvsJayGWckt83r8KpuVZ1r19oqW8eVRsAuunx0Uyx2CZSFip8xUkiQCgQCBQICOjjfBvlzIPJVKAUUXlpmZmRLJxt3GvSJgDj9moCu4uuWg6zgOKysrTE5OupmlrutMTk7uq7yjN7x+ZF5r9918dr1aDdPT0yQSCQYHBzlx4kTFJpkkSQSDQebn5wkGg5w6dQpZlkkkEiSTSebm51gILKCjo2s+dF3Dki1yZo6CalBwCqiSii5XVzSLx+PcnLhBuj1OZ2snmqqRt3LczozTEmym0+lCkiRUSUWRZBzZ4UT0JG8k3iBjpVElBU3RCMgBCpaJ6Rj0+Hvc8z/YcpBHon+PucwctmPR5T+A3/aTSCSIJWO8bLxMQcoVJ8oASZZRJQXTEQLnStFRwi42y1r9UbKFHBkzVxROlxT8SoD3NL8HAL/qo93fynp+A58KESVEWA0WFd8C6ZJrb/Y18dPd72Uxv4BfCVCwCyxll7BsizvWDNdi13j4wMOMJm+gSSrappzmurGKvoNFEsBAeICrG1eBIndZQqI31Fvymt02vqoJmV+6dAlN01hbW2N6ehrDMNzyhMiMA4FATdnwu5nu2yh2o6nrjbW1NSYmJgiFQjzwwAMlq2m9ILobWpqXESH8yCpFo11+Hcdhfn6elZUVent73Xp1pSZZoVBgcnKSVCrF0NBQyYPu9/tdcrm9bhLPx1AslXwhT9bMMrE+wVz0DqqmoWka7XonnZEDJRmVyLINw6Dv6CHmrKKHGIAu+8jaGaJaC4osE1KKXOSMlSGqRxmOHuNIUz95M49f9rOaW+F2ujiVd9A/TKvS5oplm6aJLukMhgZBAtMpMho6Ozvp7OwkvhbnWuJqkYK2OerrU3Q0U8eQDcJOCNWRsXQHSbGwbQdNVtEVDdO2kCQZv6pzJHyEdWONmcwdChh0BlvdcoQsSWTNLAvSHLIhc9p5wH2GQ2oIO1fcVazl17AdB03245dUslaaifQYaTuFLulIjoSKWjxXaefnQiwEk6lJNEnjTOsZ2n3tJa9pJNtAkiQURXH5rFB85vL5vMueWFpaIpvNoihKiXZuOBzech7VTCnfzXTf5uEFxXg8zvj4OKqqcvLkyS16nuWvryXqAWlRI71y5cq2jIjdfPZOr11ZWWFiYoLW1lYOHDhAZ2enm92WN8mmp6dZWlqiv7+fkZGRbc/xZNN7+NHGa6TlNLqmcV/4BJaUR3IUTMMkl8+ykJrjzq1pZBRCoRCGYZDJZDh69CidnZ1krDRziWlspzihZjoFNEmnN3CQjfw6S/lFLMckoAQ5GOwDQJd1d+vaHe6hO1zMbi3Lcq2YBgcH3abhhrHOaOa6O+l2MnKKiBalN9jDjcQbWGzq7lIcF46GI/gVPz1aFwvJBXyyj6SUICvlQYJWtRmf7COohDCxuBZ7ndnsDIadx3IsUnaG1kAzkiMRMxJFsHZMbucn0RI6J5qK5Pp2vZOEP8FSfhHTKeryigptVI8iydCqt7GWXy2KA22ep9+uPokpQpZkTrec5nTL6aqv2e8RYEmS8Pv9+P1+2tvfBHzTNN3yxN27d0mlUti27TpKRCLF+nylRpq3zPFOiXc86NaT6WqaRjKZ5M6dO9i2XUJtqhSqqtYlSlMLSHv9yGRZ5r3vfW9N2UW9mW4lcR8hYK7rulsvvnXrlsvN9PIeFxcXuXPnDj09PVy4cKGmH6NfCXCh7X0Ydh5VUrFxmMvdQVdU1ICK6lfQbY3uzkPElxNMTU0RjUZpbm5menqaqamporV3NEBST6CpGqqiMBwZQZVV2v3trBVW0Chq6V6Pv053oIe0lSIgB+gJHESTi9tucf5Clk+cv2HnGY3fKOryKiHyVo5ryR9xPvoguqOjSMWJMSQJBxtN1unwtxHVI6iSipyVcHw2YcLotk7OyhUFfhyHRCpBghRryXVMtYCNjS5rKLKMYRY25RttdMWHYssEpCDTmTuMRI8XjytJDASPcihwmNn0DJfWLiHLxecpoOtEtSjNWgsSEuvGGn7Zz3DwOGOBW7wev0qL1kJv4FAJa6Ge2G9ebbVQVZXm5uaSHVS5o0Q8HudHP/qRW5740Y9+xPz8PP39/XWd1/r6Oh/+8Ie5c+cOR44c4etf/3rF4aaPf/zj/N//+3/p7OzkjTfeqPv9215vXa9+B0c2m2V1dRXTNDlx4gStrds3d6Cxma4AgqmpKVdy8pVXXmnIZ5dHOUDncjnGx8ddDq3Ykgnx69nZWbLZLNFolEKhwMTEBE1NTZw9e7bu5ockSfg2+Z+SY6NKKoZdIG/nyFgpbMtmaWWZ5kKbK7spwktTUuIayUyCvGFzOzBFOBLmjm8CvxYoeqk5DvHCBgUnR1AJk3QSJM0EBznC+Ng4kUiEc+fObanfpwop8nauuP21CgSUIHk7h6M6FKwCITXkZqiq5C9OmikBdKl4HwJmCDtgIiEXmRKKD9uxUdGQgjZhC0JymI3CGhbF5l6TGnZLEGwKSToWxZo1pYu6JEnoks5AeBAHhxvxGzg49AUP06Q1o8oqQ9FhCo5BUArxo42rJH0J5AKsGEtk7SzD4fo96+DtI2AOWx0lNjY2OHPmjJsVj42N8cMf/pDnn3+e//Jf/gsPPfQQX/rSl3b83Fq0dAE+9rGP8Zu/+Zv8yq/8yq7ev13c86BrGAaTk5PEYjGampoIhUI1AS7sTsi8kpml8COLRqMlKmgCSPcj0xUSj4JD651g8zbJDh06RCKRYG1trWRUVDQXo9EooVBoV9tOWZLp8vWwkJslZSYw8xZ2DlqbWlE1ZUsm5KUp9VAsETiOQyaTIZFIYBgGuXQex3aQFRknYGITQJJldHSWYyukYzlOHD9ZUQTFcRxmszPk7SzyZkZrOHn8UgCf4sOv+tFk1fWVK1gF8k5+cwTaJpVKoqkqmhpAQaFd7UCSJYxN4M8UctxM3iRPjqgvTME2iOhhfKoP27bRZYl8waBgmkimxZq5RrvUweLCIpFIxKVa2Y5NxkpzIHCAgfDg5jCHRdzeKJYdcNAlHzkrT8ZKozs6PtmH7ujMZ2cZDA2VMBVqjberwpgIRVFQFAWfz8cTTzzBzMwM/+bf/BtGRka4c+dOTZ9Ri5YuwE/+5E9W/Mxa379dvONBt1p5oVAocOfOHVZWVjhy5AgjIyMsLi7WJae4VyHzcj+ycjGceuu0tZY6JEkiFotx8eLFEkv1Sk0yAa6JRIL3vOc9tLS0YJomyWSSZDLJ9PQ0qVQKWZYJh8NEo1G381zTYmErZJbyZKQcLeFWgu1FQM9vipZXVrQtvRaR8aSSQ8xnZ9FkHcPMk7cdnILDamKVgllA8RWzo1Qq5bItvM+HYedYzM3j4GDYRjHbdWyON51ElVXafO10B3pYyC4Uj43EmZazLMYXWMktEwoG8WkBmrRmEmYcdVMBTMfHan6FldwahpVHkcFxNBRZodnXTMYsPnOKLNMebCVuJHFMGG4ZoVlqJpGNsRRbwFwuYGclaLLAV3z+QlqYoXCxvNIst2FSAIpNtA1rvagm6bmFtdDHqsXbVWGsWghZx1Ao5Ooi7BS1aOnu5/vhHgDd8hA6AQsLCyWAA/XLO+62vCAEcfL5fMl2vtrr6/ns7UIA6OjoKJqmlbhVlDfJoEgVu3v3LkeOHCmxoFdVdQsFyLIsF4jn5uZcLqZodAgwFj80wXmenp7mwMEDtDW1ocnFrb7hGATkYN21x8HwEJqssZZfpSnQhJk1Wc4t4Q8GiAbCBOUQnWYX6VS6pDMuzi+px8lbOXTJhyM5mE4Bv+x3BwwkSeJc2wWWc0vk7TxKTmHhxl2amlvp6z2CLVtE1CiWY5FMxYvcXVnGsPIkzSRhNYxf9pN10pvOxiF02UdBLmDYhtscDGlB9EIAXVeZMafI+7JIkkRHayc5K0fBstFMjUKuwEpmmfX5daJmy5siNSEfjmITlEP4JD9xOYZhGxTsAr2BQ7vKcsV3/E5yjajGXtirlu5+xzsedAVQ2LbN/Pw8MzMzVafIdgOi9bzecRzW1taIxWIu/atRjISdygveoYpjx46xtraGqqpb+LYAy8vL7gDIhQsXavqhKYpSsdmRSqVIJBLcvXuXsbExbNtG0zQymQzRaNRt2CUKcRZysxScAn45wMFAX03XXXIPJJkjoQE66CruIDSd+/regyHn8CsBuvzdRcDxNLSFrkAymWQhPY+tOeTlfLH5VZRrwC8FSo7RqrQxfnucfD7PiRMntjBbHMeh03egaDWEhC7rWLaDrkj4VB3Jdig4+WJZAHtzsZGwHQvHgXQhhaEXWDWU4gSZrAIOsUIMn+wD2SEYLB6zYBsEW8N0S71inEdoAAAgAElEQVQkk0kmE+PMxKeLMpGoHCr0UZBMgqEQHaFODgWrWy7tFG/X8kI1YZtkMlkRdPeqpbtd7PX9cA+ArmhQ3b59m46Ojqr2OLB/tupex2FVVXnwwQcbrjRW7bXZbJbx8SJAHDt2jGg0SiKRcNW9vMMNgiYnxM73MiEExYUgGo26DJBsNsvY2BiFQoGDBw8Wp9TeeAPLsggGg0SiEcKRMNFwFE2u/9iFQoGpqSlisRjDw8M1EeO9ugJO2mQqZRc9zuw8pmXhzwR47dXXsG2bUCiEbdskk0kGBwc5cOAAkiSRMlNkzQwhNURQDSFJEoeDAxzw92I7FrrsI2VkWMot4ZN1TNtGlXU0WWM9t05EixTF2R2JRC6OZdmofhsbB1UCbAcJmYJdICSHMBzjTVEgTCJqhKAviKWZrNorNElNKJJCzsyxqN6le+0gvoKPpewKq8q6W/4RdeJagbSR9jpvBRPCtu26j1GLlu5+vh/uAdAVClfVbHq8US/o1qLRILLr3t5ezp49W5dGwl4y3UKh4HJQjx49Snt7u1ujlWWZ9fV1JiYmiEaj6LrO9PQ0tm1z/Pjxik2mvYRo2IlzKR/w8DbD1lfXmZ6aplAoFIHYU5qo9v0J4Z+ZmRn6+voYGhqq+R4X7ILbOOv1H2Itv0rSTKLKQfy+AA/0nkM/6mN1dZWxsTGCwSBtbW3Mzc0xPT1NPpplNbCMIisoisz9Lec4FOorck6VN/mx721/H2/EXmc1v0qr3kaT3sRifh5VUilYBRKFBKZlIiHR7Gtx3Sg0Wcew88UdgOJHl3XSZoY7xm1kZHp8h2iiuVhnLySRHAld0YrOFKqfpJ0k2hxleHDY/S6EYtjs7CzpdHHyzStSU03QvNFAudPvsZ7PqiZgXm/UoqUL8JGPfIQXXniB1dVVDh48yOc//3l+7dd+rer764l3POgC9Pf319TZ1zStLnPKauH1I2tvb3eza9M0922YQrzWy/M9fPiwK/Pordn6fD7Onj3LxsYGU1NTZDIZNE0jFAqxuLjoZqflAtb1hhcMDx065KqRlYe3GSaaEELbVfAwZ2ZmyOfz+P1+F4Sj0aibyTc3N3P+/Pm66oRZK8N0dtJV1mpSmzndfJaUmcABolqUQt7k9bHXsW2b+++/v6TZmTEz/Pnd59AsDcd0yOXzvJj+PkezIzRHmt1z9Pl8aLLGA61nAYgXYlxZv4hf9gMSyUxys8YbJKyHMWwDTdYYCA0xm72DJuv4JZWeQA9xI07OyRHRmrBti7XCMrngEUKEUSSFznA7CgpIsJreQEMrqeFWq8en02lSqVSJoLl3+CASiTS8vPBWuEbU+/zWqqX7R3/0R3W9v564J0C31qhXk6BSCD+yUCi0xZdtN8Lk9WS66XSal19+2eX5btckW1lZYX5+nr6+PtfbKZ/Pk0gkSCQSzM/Pk8vl8Pl8JSC3nVSfN0Qm3dLSUpEPu1MIdkEwGHStZ4SKVTKZZH19ndHRUUzTJBKJuNKSkUik5nn9+dwMkgM+ueguEjNjRLUWmvVWbNtm+k5x2k7sFMojb+dQZAWfWszYghQ5vX09fThZZ8t9FOCV92cBiULBLAoD+XUiWpSQGiFlJgmpIU41nSGiRegO9G5yglUkSeJ76b9Cl4tDGigKWStLwo7R5G9GlqCVFlaMNWRHpjXQRH7Jwt/uf9OxQy6qpokmHxSfM28ZSNzrTCZTsuglk0lisRgbGxvutdT6PJTHfpcXUqnUO1JLF37MQHc3D4/IItPpdEU/sr18fq2gG4vFGB0dJZ/P8+CDD7pAX6lJtrKyUlLf9j74YgTTW/wXIJdIJFhYWCCXy6HruvsjLQe5TCbj+sKdPHmyZk+4WkKSJHRdJ5VKEYvFuO+++2hvb3cXi2Qyyd27d8lms+i67gJDNBrdQg8DKDgGmqS5n43jYDmmK5spGonVsruQWjSfLNgFNFnbVOpSaQm1oEbUkvuYz+fd+7iytEwiGEeyZfy6H0uyCUkhHmx9/5ZjyZt+ayI0WSVn5UqzV1nftPuBVl8HYa2JWCKGWTAY7D9MS1NLiRqcVyTePc4mGIvje3cfIq5fv05XVxeO45RIN3qdJSKRCMFgcMeMeL9NKePx+LbTpG/nuCdAt1GandU++9q1a3vyQ6sWO3FvM5kMY2NjWJbFsWPHGB8fx+/3V1QAEyaZQhJyO2dkbwgg9s6wVwI5VVWxbZtCocDAwADd3d0Nve+iZHP79u0to8c+n4+Ojo6SczQM402QW1khk8mU0MOi0SgBOUTGSqGjF+UaLZvbt6bQHL2me6TLOhfa3sfltZfJ2Tk0SeNC2/tQ5a0/G5/P546ZO1k4fuAk88xgmRaSJRNcjXBp6tKOfOfh8HGuxq6QtbIAhJQwB3wHkKQi1SKTy7C+tk4wHKS1vYVWrWXzb5R8jheABS/btm0XiB3Hca2T4E2pUbH7KJduFHXi6elpMpmMC9pi4Su/jv0uL7xTZR3hHgHdeqMWIXPRqEqlUvT29nLw4MGGg3u1TFeoem1sbDA0NER7e7tbzy1nJOTzeVeha3h4uCFbLi/ICTWy6elp2tvb0TSNlZUVpqen0TStpDRRKdusJYRweiAQqHn0WNd12traSpp2XnrY1NQUqWwK2kyUoIxt2dirMHzweMVSQrVo93Xwd7r/Lnk7j0/2VeUWCzt6b/Y8aA9hOgV8ss8Fxmp851Ao5N7LM00XSFgxVFmly3cAVdYwDIP1mTiELdq72op6BXKb+7nl4S0tiNjOOsk0TfL5fMm/efntggEiwrKsLSI1ArTD4TDZbHbPpTwRpmlueSbeqa4RcI+A7m7kHavVIL3DFYcPH6ajo4Pm5ua6jrFbc0rbtpmZmWF+fr5kYEFkKkICUYzm3r17l9XVVQYHB132QiNDyF62tbWV1JBFeLPN5eVlt2HnZSSEQqGq5yVGtE2nwNDxAYLBEAq7pyyVg4PjOCzcXWDq9m3C4QiqqjI5Ocnk5GRJaWKn6TpZkgkolTVrc7kcY2NjOI7D6dOnXdF5KIqTq2U/se34zkLuMJlMuo0uI2JiGAZra2vF77mjDVuyUVDrHi6pZp2USqUYHR0lFArh8/ncZ1L8f1EjLq8TNzU1lWSbtm2TyWRIpVIYhuGKKZVr6NZbJzZNc0sZ652qpQv3COjWE9VA1+tH1t3d7Q5XJBKJuhgJollX68CBaIAtLS1x+/Zturq6SgY7vE2y8+fPE4vFmJubIxaLoaoq4XCYRCIB4DY+9hpCUF2WZU6dOlUCJN6olm2K0kS1bb/f7+fu3bvMzc0xONJPpDkMkkQBAwsbHd+O48G1XMOtW7fw+XxcOP/ekkxJZGmiESakBL3b/kgksu32WLBI7t69W7URV2uU853hzenCiYkJZFlG0zRu377N3bt3t9DsdrvY2rbN9PQ0y8vLJW4h4m/V6sTVGnaidBIOh5mdneXMmTMArsWP14FY07QS5sR2u6RK9eF3QfcdFOVc3Up+ZN4f6G6EzGvt3CqKQiaT4dKlS4TD4RKucaW6raiptba2curUKdeWPZlMEo/HmZ2dLaFdealhtYQoqcTj8V3XrzVNq7rtTyQSjI6OEo/H0TSNzs4OZE2iYBTQdR9IEjYmDvquQdc0TW7fvk0sFuPYsWMV637VsrR0Ok0ikWBxcZHx8XF3YMILcpqmuQyW9vZ2zp8/33CNAaEDHIvFOHnypAvGYreTSCRKvm/BnKjU+KwWsViMW7du0dHRUSJ9KaJaeaLWhp13t1fJ4scwDBeIV1dXyWazJaDtFQC6l1wj4B4B3Xo1dQVXV4zO+ny+Ej8yb+x2dHinuqRgQySTSc6fP++yIaqB7fj4OLqub8k8xQMtOunVfpjbAbEY8ijn/jYqNE0jEAgwOzuLpmk89NBDm02nBAU7z/r6BoZRFKAJhoNIRpxIqD51M6+Gbr0DFFCqcNbbW7SxEdtl0awTHnKyLNPd3U1zc3PD1bTEcQ4ePMi5c+dKrsHrUSZodlDa+BSiTqLM4623y7KMaZrudVQT79/uHsHODbuNjQ236VqpYQfFXVJ5ndg0TVfac35+nnQ6jeM4GIaBz+fDNE3XCj6RSHDkyJFd3eO/6ZB2mOzYvWTRWxjiC64lxsfH8fl8rK+vY1nWjs2n6elpFEXh4MGDNX3+66+/Tn9/f9XPFHXMeDzO4cOHWVhY4OzZsxXBVrw2k8kwNDS0626tF4jFj1MAsaqqxONxOjo6OHr06L5kbVNTU6ytrTE0NLRFVtMgh0VxUsuybQp5g8RKMeMsVzcT9ddyIE4mk9y6dYtwOMzg4GDDRllFOI7D7Ows8/PzDAwMuEpm4l4ahkEgECjJNust8+RyOW7duoUsywwPD+95msu7u0gmk6TTaSzLwjAMOjo6OHToEJFIpGHDEFAEzfHxcTKZDCMjI/h8vhIw9kalOnGlsG2bK1eu0N3d7daLn3jiCdLpNKdPn+aXfumXOHPmjLtQbhd7FTB/+umn+W//7b+5GfsXvvAFfu7nfq7a4aqu+PcE6IrVcKfI5/O89tprGIbBiRMnqvqReWN+fp5CoVDzqnrjxg26u7u3fJneBl1/fz/d3d2Ypslrr73mgq7Yktm27cpSDgwMuDq4jQwBVI7jEA6HyWQyrnHgbkoT5VHu3tDb21vxx1W0nSlgYyEjo5aVFkS3X4CH6PaL7WcsFsMwjIpmno2IWCzG2NgYra2t9Pf3V1yUvIuaOFexqJUDcfn36K0Nb+eRt5fI5XKMjo66GbqYBCxXiqulll0tRIZ++PDhqnTCauUJEdWA+PLlyyVZv2VZ/It/8S84ffq0u/h9+ctf3vEcP/vZz9La2uoKkG9sbFTUwv3e975HOBzmV37lV7aAbjgc5tOf/nQtt6TqD/aeKC/sFF4xbyFkXuvDraoq2Wy25mOV14C9tu7eBh0Ut1qZTIb5+XmampoIBoOuHGJvb2/NNjn1hGEY3L59m2QyucVk0gsesViMmZmZLVlcLUAsyjbV3Bu8ISGhUb0UU6nbL77P6elpAoEAtm1z48aNLRnxXrJ2wzCYmJggl8tVVBoruYYK235hwiiAWEyuicETMWV3586dfasNiwx9YWGhKqB7leK8texgMFjSVKxWLhMsBcdxOHPmzLbPxk7lCe9Oz8uc8CYk4v25XI5f/MVfZGhoqOb7sVcB80bFPQG61bJAr06B0NZdXl52RUBqiXqFzL2vF9zNpqamkgad9yE7efIk8Xi8pMHU0dGBpmlks9ldc1/Lw7Zt5ubmttDRvFENPKoBsZej6/P5XM5wLpfbF2EdeDPzbGlp4aGHHnKzskqMBJHFV9L7rRaClzw7O8vAwIBr2llveE0YyyfXNjY2mJ6eJpvNuo050zRrotnVGslkktHRUVezotp1V2JOiFq2aHJNTU1RKBRKvvNIJML6+jrT09MMDg7uSuJQHB8qA7Hgqwsam/d3uLS0VHcjrREC5L/7u7/L1772Nc6dO8eXvvSlXTWb7wnQLY9KfmTix7lf8o4iFEUhnU7z6quvIkkS73nPe9wsqVLdVtd1NjY20HWd973vfW6ToJz76t3y16o9II65urrK5OQknZ2dNevniqgViAXtqr29nb6+vj3LRpaHAPRqGrfVGAkCiBcWFnakhiUSCW7duuUuko0W4XYch/X1dbemKKQjt6PZiQWj1qail/mw25KLl0VQSaBodXWV69ev4ziOKyXqOE5duhg7HV8kLAcPHuT48eNuTTiXy/Gf//N/ZnZ2tmJWvZ8C5p/4xCd46qmnkCSJp556ik996lN89atfrftz7hnQFduQan5kIvZTyNwwDJaXl8lms5w6dcpdBSuBrVjFU6nUlm1+OeXKMAwXiBcXF13tAS8QV6oXCgM/4f7bCA6vOH8BxLIss7a2xqFDh+js7CSdTlfMiMX/1QvG3gx9cHCwrvp2tSyunBrmnfLr7++ns7Oz4YCbTqcZHR0lGAxy5uwZbMUi52Tx4a9IsxMSjYlEombLJPHs9/T0bGE+7DXEd766uko8Hnefb6E54R0Zr2dApjxEMy6bzZYMm0iSxNWrV3nsscf4+Z//eaampiqWrfZbwFzEr//6r/PII4/U9X4R9wzoiiylmh+ZCFVV65J3rAWkhYj54uIira2ttLW10dLSUjLz7m2STU9Ps7i4SH9/PyMjIzs+kLqu097eXkLA9yqGCaEaoRgWDAZZX18nm81uaxe0lxCUN03TSgA9HA6XZMTZbJZEIuFuqesB4vX1dcbHx135zEbUPEvML3t6uHv3rms1r+u6W56wLKsiR7feEM/G2toax44dIxQNsWQsYJsmDqBJGl1aD3KZxU49lknBYNDNjr07q0ZGKpXi5s2bW8oVlbQ7xKSiN3P33vdqLJS1tTXGxsbo6+sr+V3k83n+03/6T7zwwgt85Stf4dSpU7u6hr0KkAvABvjGN77ByZMnd3Ue9wR7AeDmzZu0tbXtCDD5fJ5r165x7ty5mj7XNE2uXLnCe9/73i1/806x9fT0cPjwYTcTOHr06Ba5RdHN7+7u5tChQw1tnAiAEw1DkeE3io0gwjtAUat7Q6Xz9NLXDMMoETT3+XxMTU3hOA7Dw8NVJ+L2EjvRzLwcXXGeYjTXW8veDohF5im+b1mWWSssk7ZTrvpZwTZoUltoUmtzqPaGqD/fuXOHtrY2Vx1MlFC8C8ZuM3fBpFldXWVkZGTXyl5es9Ny5kQwGCQWi+E4Dvfdd1/Jjuy1117jscce4x/+w3/IZz7zmT3RAdfW1vjlX/5lVwz/j//4j2ltbd1WwLyrq8sVMP+n//SfcvXqVSRJ4siRIzzzzDMuCFeIe5syBkUwqEVgw7IsLl++zIMPPljT5zqOw0svvcRDDz1U8u9iVW5paWFgYMDN1tbW1piZmWF4eBhd1xHOvBMTE0QikZLXNirEyKgYI+7r60NRlC38XGFjvpstf7l7Q09PT8O2rwKI4/E4c3NzJJPJkjq2AI9G3DfTNJmcnCSRSNRd8xQSn2LbL0bEBRCLc3Ucx/WLO3bsWAmILBnzFJyCK91oOiZBOUSbVt9WV5QrQqEQR48eLQFVUULxUu3KF4ztGAkiRIO3s7OTw4cPN5xJY1mWu2iEQiG3zDM5Ocno6ChLS0uMjY3x1a9+dddZ5d9g3Puga5pmTSyDaiC6Xbz44ovu61OpFLdu3UJRFIaHh90yhigjGIbBnTt3XK6mkLg7cuSIy0poZAiFLr/fz9GjR3fMZCtlmpZl0dzSTCgUIuD3E4mUZnCCMdDc3MzAwMC+uLyurKwwOTnJgQMH6Osr2uF4zzORSLgWP16AqxWIyyfWGrVoeK2IEokEq6ur5HI5IpEI7e3tWxa2uLlB3FxH3cx0TadAq9pJWK0N/L0c7nK9hJ3OsxyIq1kmWZbFxMQEyWSS48eP70u5olAouII4YpBCXN+3v/1tvvzlL7sjxbFYjN/7vd/j/e9/f8PPYx/jXdD1hhdEa3392bNnGR8fJ51OlzzslZpk3vn/Q4cOIUmS+6MUtUJvprmbMkM+n3en1YaHh/cm6Cw7mGYBwzDIZnNMT824o5f5fB5JkhgZGdmX2nAmk3Fr8UNDQ9s2+7xuB/UAsVgog8EgR48ebfjCB7i6Ei0tLfT395c0P0um1qIRtFYFx2ciyTJRuZkmtaWmBUDoJTQq8/QuwOKeZjIZCoUCLS0t9Pb2NsTWqTyE1VV/fz9dXV3uZ+dyOX7nd36Hl156iWeeeYYTJ04Ab0pR7rele4PjXdD1Rj2ga5om3//+9/H5fAwMDLgPSaUmmZeMXm0yx9s9Fw+7bdslwwfbjWcK+UfRiNstj9QNmS2Ph21aTN0uNgY7Ojpch9zyBWMvtULvePDw8PCuxeG9maYADgHE4XCYdDpNJpPh+PHj+7JoiHJFMplkZGSkKjfZS7nylnq8Nfdq48OFQoGJiQmy2SwjIyMNdevwHmNsbAzDMEoWjWQyWWJHtB1bZqfwDlKMjIyULI6XL1/mk5/8JB/+8If55Cc/+U4D2Epx74OuZVk1U7tefPFF3ve+92370HjFuwuFAj/xEz+xrSfZ8vKyyws+fPhwXdmrADXxY0ylUkiSRDgcpqmpyWUkCKK62II3pMbmBV3HIZNJs7BwF03R3eaP9zzLa5pe3mstAwhehwghDt/oWqGgmYlaoSDaezPinZpgO4X3OnZbrthpfDgSiWAYBrOzs/T397u83kaG9zrKM09veKfrEolEzZZJIoR0afkgRTab5Qtf+AKXL1/mmWee4fjx4w29vr/BeBd0vXHp0iUeeOCBqj+61dVVxsfHaW1tZWBggKtXr3L69Gk0TdtSSkgkEoyPjxMMBhkcHGyY7bRXc2B1dZVYLIaiKLS3t9Pa2ronp4aSkAAZ8rkcq2urBAJBIqEImlobIAkgjsfj7oIhJsHKR3IFb9jn8zE0NNTwhiKU6uh6j1Fee/XWNOsF4mw2y+joKLquN/w6xPjw6uoq09PTbk+gUka81+9eiOyI/kS91+EVsU8mky5trZzdMTo6iqIoHDt2rOT+Xrx4kU9/+tN89KMf5bHHHrsXsltv3PugW4/S2GuvvcaxY8e2bNMEjUjTtBKq0quvvlqi/CRJkuviUCgUGB4e3peRV+9Y7fDwMH6/vyTLTKfTJV3+pqamun+MhmGwcHce3V8UJPdp/j1/6+WZu1A1g6LVdVdXV0We5l5ClCvW19drprJVA2IvP9cLxF7R772URHY6J1GiGh4eprW1dYuOQyKR2NOWX+zi5ubm9izAXh5edbPl5WWSySR+v99NFObn5xkaGuKLX/wir776Ks8++yzHjh1r2PHfRvEu6Hrj2rVrHD582G0+5XI5dwrGK3wtarZTU1MsLi4SDocJhUKk02nS6fS+qUIJRbKlpaUdVca8DRux7RNDEt4fY3mUazHs19bVqzQWiURKMiNBmBfnWY9+rjcE86Gnp2fP5YpyIBa0ME3TyGQytLW11cQS2U0kk0lu3ry5raKZN7wDMslksmRSUdzX8rHcTCbDzZs3XX7yfmSX+XyemzdvusmL0ISOxWI89thjXL9+HU3TePjhh3n44Yf58Ic/3PBzeBvEvQ+6tco7QnGQoquri2g06g4THD161AW3Sk0ykUndvXuXYDDoNu28oLHX7E3Y9kxNTdHT07OlplprlHNz8/l8CTdXLCRtbW01/bh3E7Vo3HpLKF79XO+Csd0IaTabLWE+7AcQGobB2NgY2WyWzs5O9956ea97mViD4n0QGst7FQryWsGLLb8YyxXlgOPHj+9bli5U8oShqohMJsO///f/ntdff51nn32WQ4cOce3aNeLxOA8//HDDz+VtEO+CrjfGxsYwTZP19XX6+vpKsqNKTTIhGNPR0cHhw4e3KFt565kCNLwNsFoySCGHKECq0XXCXC7H6uoqMzMzFAoFdF0nFAq557kX0PCG6LQLal29gitezQFRQhF1Qu9U3czMDCsrKxWF0RsR3mGQSmpj1QYlytkdO91TYRW1X47TUBynvnnzJj6fz5Uq9eojNKI/kMvl3GMMDw+XNJ1feuklPvOZz/Crv/qr/Mt/+S/3ZZF/G8a7oCtes7Kywo0bN4hEIpw+fbrk4ShvkiWTScbHx/H7/QwODtYkGGOapvsjjMfjZDKZEnGapqamEt6jqA0LC/X9qA2bpunO/wuQqrSN9lLCmpqaapJCFOGVRGx0ucJbJ1xZWSGRSLg28WLRaIS6lQjhjhuJROraggsg9tZeq2XE+Xy+6tRao8KrOFaeQQtls+2aYLUI1Xjrw6IGLSKdTvP5z3+eGzdu8Oyzz3L06NGGX+PbOO590AXcZk2lEII4ovkgyzKHDx+uCLZeGcGhoaG9DR5Q3KKKbNjbBLFtm1wux+Dg4FtSU63m3iCiEoe4nIlQiUPs1bjt7+/flzqh6LRLksTw8DCyLJdkmeWL2264pF6QOnbs2J6/d6gMxLlcDtM06ejooKenp2G7DG9sbGxw69Ytt0xVy33wLm7innqdnIVimPj+s9ksN2/edIdOvAnMD37wA37rt36LRx99lE984hM/LtmtN348QNcwjC1eTNlslvHxcfL5vPtDWlxcJJVK0d/fXwK2QhFqdXWVwcFB2tvb9wUIhZNEc3MzqqqWdM29oLEX8PK6N+zFN6wah1jYZm9sbLijnPsxLioYA0tLSzs2LmtpKlabrvIaQu7XNt+rl9Dd3V0CxuUZ8W6/f6804vHjx/csFiSAWDwDwphTMHgGBgbo7u52gTiVSvHbv/3bjI2N8eyzzzI4OLin44v4zne+w2OPPYZlWTz66KM8/vjjJX//H//jf7guEOFwmD/4gz/g9OnTNb13n+LHD3TFKO7q6qpb1BdNMjFkcPjwYZqamlBV1d0a15IR7jZisRjj4+MVgdCbEXkHD2qdVBPhpZkdO3ZsX8oVQgt4aWnJbSqW110bwSEWQtZdXV27HnsVFvXeXYaX8yoUzRplCFkpRONSyDtWmowr//69E4DeZ2A7IK7Fp2yvkclkuHHjhrurSKVSpFIpnn32WeLxONevX+cf/+N/zG//9m83rGQiDGT//M//nIMHD3L+/Hn+6I/+iPvuu899zYsvvug2CJ977jmefvppLl68WNN79yl+PDzSRLY6NzfH7OwsfX19PPjgg1uaZNFolJ6eHlcpLJvNEggE6O7ubsiWsjwEJc00zardaTGBFg6H6enpAd50PhDqW94sU9QyRd3NOx68X2aWUKpx+xM/8RPuttG7NZ2YmChxvBDnWut2X9Q7LcsqEbLeTZTrvQrOazweZ2Zmhng8jq7rhMNh5ufn3XvbqEam2OYfOHCAc+fOVV04qn3/ou4u9Aq8dXcBxrZtMzo6CrCjT9luw8sfHhkZKeFBJ5NJIpEIGxsbfPzjH2d6epqf+qmf4q//+q8bAryXLl3i6NGjDAwMAPCP/tE/4k//9E9LgE0j42IAACAASURBVNM71v/ggw8yNzdX83vf6rinQHd5eZmxsTE6OjpKLHrK67aiTnX37l23oWbbNvF4nPn5eZLJJJIk1UxdqhbeBpYoV9QTlZwPLMtys6Hbt2+TyWRcdbPW1lZOnjzZEI+t8sjlcoyNjeE4DqdOndoChJqm0draWtJI8W73y4XWK3GIvc64wiWi0SG2xdPT07S1tXH//fcjy3JFH7jdKppBcREaHx8nl8vteuHw2uZUAuKlpSVu3rxJLpcjGo3S2dnp1mEbWVdPp9PcvHnTtTESC63jOPz1X/81TzzxBJ/4xCd45pln9mWHOD8/z6FDh9z/PnjwIBcvXqz6+q985Sv87M/+7K7e+1bEPQW6wpFU/JArNckMw3DVuYaGhkq2el56kwC3eDzuglv59Fe1+qCXryi2NI16GBVFcR0FhHuDoih0dHSQzWZdYRQBbiLL3G32U09NtTwqOV54OcRzc3Muh1jTNGKxGB0dHfvijAulQFjus1bJB04ocK2trbnmjDtRwrxc6/0YOhFArCgKS0tLtLS0MDg46HJwl5aWXEffvYoTOY7jfvflKnPJZJLPfe5zzMzM8Gd/9mccPny4YddY6TzKo9o9/au/+iu+8pWv8IMf/KDu975VcU+B7oEDB1yObTnYCh3S5eVl+vv7OX78+LY33wtuIrwshPn5eRcwBLCJGtf4+DhNTU072o/vNmpxb9gpc6tFZ8CrcdsoO/hyh1wxvZRMJmlpaSGRSHDp0qWGNhW9LI5agVCSJILBIMFgkAMHDrifI7JMcW+8DTC/38/8/Dx+v3/fvvtKY8KAy8rxZsSiRry0tMTExERdKnHCnqelpaUkaXAchxdeeIF/+2//Lb/5m7+5b9mtNw4ePMjs7Kz733Nzc+51euP111/n0Ucf5bnnnnOTg1rf+1bGPdVIKxQKmKZZMkkGuFnnXqa8KoUYOojH46ytrbGysoLjOLS0tNDa2kpTU9MW88C9Hm+37g3VeLmVFMLq0bjdy7XMzc0xNze3xXBSnKtY4ERTqR41MxFeQ8j90NJ1HIdkMsmdO3dYX19H1/USA8l6znWn8PqUDQwM1P2ZXiAW9fdya59QKMT8/DwrKyscP368ZPeXSCT43Oc+x/z8PM888wx9fX17vqZawjRNhoeH+e53v0tvby/nz5/nf/7P/+nq7QLMzMzwwQ9+kK997Wsl9d1a3rtP8ePBXvjMZz5DOBzm3LlznD17lu9///vouk5fXx8DAwP7knmYpukKrRw9etTd9nsBo1rzq57YD/eGcoUwMTIM0N3dTXd39671ELaLeDzOrVu3XKujWsDDCxjxeLyimpmX3VFuCLkfWrrwpnh5W1sbR44cQVGUHTWThRxiraDpZT/sxaes2meLc11bW2N1dRVFUdykIZlM0tPTw6uvvsqTTz7JY489xq/+6q829JnYidI1OjrKP/gH/4Bbt27R2trKJz/5SZ588kn+63/9rzz55JP09PSwsLBAPB53bX1UVeWVV14B4Nvf/jb/+l//ayzL4uMf/zhPPvlkw859m/jxAN1bt27x8ssv853vfIf/9//+H11dXZw9e5Zz585x/vx5Tpw40TDg9Wadhw4dore3tyqQCo2BeDzuTqmpquqC8Hb1YS/zYXh4eF+4sF5N1Z6eHqLRqJsJpVKphtHBDMMoEePe67WIMWwvuMmyjKZppFIpurq69s1eqFbxchGCieI9V2DHwZP99ikT5yZA/b777iMYDLpA/Oyzz/KNb3yDtbU1PvjBD/JTP/VT/LN/9s8atojVQulaXl5menqab37zm7S0tPDpT3/a/duRI0d45ZVXGqqU1qD48aCMHTt2jEOHDvG//tf/4k//9E85d+4cV69e5eWXX+b3fu/3uH79OsFg0AXic+fO7epBXl9fZ2JigpaWlppqd4qi0NzcXFJ7FZ39eDzudvaFKE1TUxOhUIiFhQWWl5cbLr/nDWFlEwgEOHv2rNuh99ayvSOjy8vLJdNfXjpYtfCOCG8nlF1vKIpCU1OTCwD5fJ7R0VEMw6Cnp4dMJsMrr7zScA6x4MMeOnTIVdHaKaoxUQQQeymBYrsfi8XI5/MuI2U/IpFIcPPmTTo7O0sobeFwmJdeeonnn3+ez33uc3zkIx/h5s2bXLlypaGNqFooXZ2dnXR2dvKtb32rYcf9m4x7CnQBgsEgzz33nPvfDz74oOv86zgO6+vrXL58mZdffpmvf/3rJQwDAcYtLZU9qzKZDOPj4wCcPHlyT7Yp5Z19b314ZmaG9fV1NE2jpaWFTCZDLBZrWG0Q3hweEc247TIXTdNoa2srYS4IWUHBIS5XMmtqakLTNHf8WtCN9iPrFPXh+fn5ilSz7RaNekaG8/k8t27dAhrDhy1fNACXZz41NYXf78e2ba5fv94QCUxv2LbN7du32djY2ALqsViMJ554gvX1dZ577jl6e3sBuP/++7n//vv3dNzy2CulS5IkfuZnfgZJkvjn//yf8xu/8RsNPb/9iHsOdLcLSZJoa2vjQx/6EB/60IeAN91VL168yAsvvMAXv/hFV/5OgHBfXx9/8id/wvnz5zl69Oi+qFoJQ8v5+XkCgQAf+MAHUFXV3eYtLCy4W1IvsNVbH/bS2fr6+hgaGtpV5iIEZ7xDB2LRWF9fd2l2UBQu369hDQHqzc3NValmlRaNSjuNajQ7b6YuJED3I4RDrmmaXLhwwd09eJXX7ty54yqvlQNxrfdXlCzEwIa3gfn888/z9NNP8+lPf5p/8k/+yb4zE/ZK6frhD39IT08Py8vLPPzww4yMjPCTP/mTjTzFhsePFehWClmWGRgYYGBggI985CNA8eG/du0aL730Ek8++SQ3b97k1KlTTE5OcubMGc6dO8fQ0FDDsk7BHU6n01ucfUXTRWQb3vrw1NQU6XQaVVVrmvwSABWJRBpOaZIkiUAg4GZnwnBSCJcvLi66/FEBFoLdsZsftremuhsN2ko7DW/2Pjs7Sz6fR9M0stkskUiE+++/f89aBpXCW1OvJCOpquoW+qJ3AlAscDuVUYRubyKR2JLdbmxs8Pjjj5NIJPjOd77zltGq9krpEq/t7OzkF37hF7h06dLbHnTvqUZao+N//+//zY9+9CO3m/rKK69w8eJFLl26xMTEBJ2dnSX14XprlV73hr3UOr2TX/F43NUX8MoezszMuHbt9Wrc1hq1CJdXsvKpR7i8EYaQtYQQrV9dXaWrq8stUTRamCiXyzE6Ouq6LOxlIfSWUYRKmBjoUVWVxcVFent76evrK8lun3vuOT7/+c/zW7/1W3z0ox9tWHZbCyvhYx/7GBcvXuSJJ57g6aefdilds7OzW9779NNPEw6H3UZaOp12F/F0Os3DDz/Mv/t3/87dxf4Nx48He+GtDMFeuHjxogvEQlxHUNbOnDlTtWnj1TAQVKNGnpvY6s/PzxOLxdB13a0f1quVu1OIrDORSOxKEtGrQSzAwsvuENm7cIrYD0NIbwihne7u7i287mrCRPXycr06tPtl+wRvumukUikCgQCGYaDrOj/84Q/x+Xx873vfA+D3f//36e7ubthx62ElfPGLX+Qv//IviUajfPzjH+fxxx/nwIEDPPbYY3z2s5/l/vvvZ319nUwm43Kgb9y4werqKr/wC78AFJ+hj370o28VHayWeBd034qwLIubN29y8eJFLl++zKuvvoplWZw6dcrNhiVJ4sUXX+S9731vifllo6Nc41ZRlC2cXKCEPxwOh+uuD4uR10Znnd6aqxfcOjo66OrqaqgojYhCocDY2BiGYTAyMlLzd1POd/bSwbz3VoD3W+FTBm8uHkI5zyuc/+Uvf5lvfetbmKaJaZr09vby9a9/vWG7oJdeeomnn36a559/HoDf+Z3fAeCJJ57Y8tryDLae976N48eDMvY3HYqicPLkSU6ePMmv/dqvufP7V65c4fvf/z4f+9jHWFtb49SpU9y5c4fz589z7ty5hkpJ5vN5xsfHKRQKW/QFhHhKeX3Y26AR9WFRc61WHxaTXoFAYF9GXkXNVVEUVldXOXToEAcOHCCZTG4ZbfZmxLsBMO+Y8G7KPMJgs1y7Q9DBZmdn3SEZIZJ/9OjRfRGuh1JN3XKxnbW1NT7zmc9gmib/5//8H7q6utxdWyNlQPfCSng7itQ0Mt4F3X0MMb//gQ98gBdeeIF/9a/+FY8++iixWMwtS3zta19zHXlFNnzmzBmamprqrg8Lha5apR0r8YdFXTAej7O4uEg2m3X1ZwVbYm5ujo2NjX2d9DIMg/HxcQzDKFE1CwaDJaI0XunDiYmJkq1+LWPYwv2g0XoJ5XSwZDLJjRs3CAaDdHZ2srq6yvT0dMM5xEKutK+vj5GRkZLa7Z/92Z/xhS98gSeffJIPf/jD7t8kSXIX4kbFXlgJb0eRmkbGu6D7FsVTTz3l/u+Ojg4eeeQRHnnkEaAImBMTE1y8eJHnn3+eL3zhC2QyGU6cOOEC8cmTJ6vyQkV9uBEKXeX0KtHVj8VizM3Nsb6+jqqqNDc3E4vFsG2baDS6L/oSlTr53pAkiVAo5DoxQOnkl1ems9y12XEcV3/42LFj++KOC6UWQCdOnNiSTTaKQywU1PL5PA888EDJsMrq6iqf+tSnkCSJ7373u67Y0H7GXlgJb0eRmkbGuzXdt2kYhsHVq1fdJt0bb7yB3+/ngQcecIHYtm1efPFFHnjggX2tDwsBHNHA0jStRAdBeKl5qWC7IfDv1hByp/DS7LwaE+FwmIMHD9Lc3NxQY0sRu/EpgzcHT8o99bxlFO8CvLq6yvj4+BYFNcdx+OY3v8l//I//kaeeeopf+qVf+v/bO/uoqM5rjT+HD0n9AJxZGSFAQRxAjAiGASUXlYpIgwZszbI2phqrRo2ItzGyvMvVaBOD1MRqom01kfjRtLHQVmgUjMYWP5LyZbAW4RIuSkBmAAEVUISR2fcPct6cYYZhYM4gA+e3FgvOzHnPnHc57nlnv89+9qCtGPtjNNMzp/sYTWrERNpIs3WICHfv3kVhYSEuX76M9PR0NDc3Q6VSITg4mFXUidnXTWga05uFpPDcnu3o+a/O/Nfs3lZswtXg5MmTrSZp41UWbW1tUCqVTDVx7949g35qvB/GQBD69orRp6ynhrilpQWdnZ144okn8PDhQ9jb22PKlCl6q+iGhgZs3rwZjo6O2L9/v6gFHX1JwYgImzZtQkZGBpqbm/Hkk09i/fr12LZtG+RyOevkQURoampCS0uLnirB2dn5cZnUiIkUdIcT/MogMTER9fX1yMvLQ0FBAQoLC1kn29DQUISFhbGNlP4GYt4v9qmnnoKnp+eANvqEX52FgU1o9HPv3j1UVlbCw8PDag0hge/8EkypLPg2PsLVcE+/5L5yvnxbHWv2KQPAPHLlcjk4jkNrayva29vxm9/8Bk8++SS+/PJLvPnmm1i2bJmo92COFCw7Oxv79+9HdnY28vPzsWnTJrYRNoQNasTGdoNuc3MzfvKTnzAD6vT0dIP8W01NDZYvX466ujrY2dnhlVdewaZNmwB0B6gPP/yQfdKnpKQgLi5u0OcxWDx69AjXr19HXl4eCgsLcfXqVXAch5CQEFbIERAQ0GsOltd12tvbW6VRI2+u3tTUhPr6euY/7OrqyoKbmJpl3gRnII0nhd0j+GAsNAIX6p2FvgyTJ0+2moa4s7MT5eXlICKD11Gr1diyZQva2trg7e2NsrIyZgAlFubIudauXYuoqChW4RkQEIDc3Fy4u7tLQRc2sJGWmpqK6OhobN26FampqUhNTWWtlnkcHBywZ88ePPPMM2htbUVoaChiYmLYp+8vfvELPTu44YyDgwOCg4MRHByMtWvXgojQ1taGK1euID8/H6mpqfj6668hl8sRGhqK0NBQhIeHw8XFBadOnYK3tzf8/Pys4i8BdEvB2tvbce/ePQQFBUEmkzHDcr7VjBj5YaEJzkBd2ox1jxD6z2o0GpSXl0Or1UKr1cLNzQ0eHh5W093W19fjxo0bmDRpkt5mmE6nw1//+le8++67ePPNN7Fo0SK9vK6YmCPnMnZObW0tW/nbmkGN2Az5oJuVlYXc3FwAwIoVKxAVFWUQdHnDbaBb7B8YGIja2trH2vFzqMDv3EdFRSEqKgrAd7rUgoIC5OXlYffu3aiuroZKpUJkZCTbAe9vsURf8AUbcrlcT2XBKxCErWZ4/fA333xj4Onr4uJiMmXCd1jo2UhRDISaXJlMhrKyMowbNw4TJkzA/fv3UV1dze5XWHhiiRSss7MTZWVlsLe3N5C11dXV4bXXXsO4ceOQm5trUNkmdnrDHDmXqXNs0aBGbIZ80K2vr2cB1d3dHQ0NDSbPr6qqQnFxMWbMmMEeO3DgAI4fPw6VSoU9e/ZYTR5kK3AcB3d3dyQkJECr1aKyshKZmZl48OAB8vLykJWVhR07dqCzsxNBQUEsPzxlypQB6VhNNYQ0hp2dHdt841dMxqRVPfPD9vb2rItHz1YzYiLsUyaUmwkDnrC0mW9JL5SCmTKuF74OX7TR091Mp9MhPT0de/fuxc6dOxEfHz8oygRz5FymzrFFgxqxGRI53Xnz5qGurs7g8bfffhsrVqzA3bt32WPjx4/HnTt3jF6nra0Nc+bMwbZt2/DjH/8YQHfQ5nf0X3/9dfz5z3/G6NGjB5QfNie/bGsIe8n15OHDhyguLmb54evXr2Ps2LF6Jj/f//73e/3qP5CGkP1B2Fm4sbER9+/fx5gxY+Dm5sbyrWJ/1bekT5mwsSkvBRMaEzk7O7McLd+w05gRTl1dHTZt2gSZTIa9e/daLRVkDHPkXKdPn8aBAwfYRlpSUhIKCgqGukGN2NjuRpowCa/RaBAVFcU2LIRotVosXLgQsbGxeO2114xea+3atcjMzER9fT1SU1Nx584dg1SFRqOBRqPRyw9nZmZiypQpSE5OhkwmY/llY+OHM7zEhzeBLywsZE0yeZOf0NBQjB8/HiUlJXjw4AGcnZ2t0hCSh69c02q18Pf3BxEZ+DX09JcYiBJD2NJGrFU0b0wkVHhotVrY2dnh4cOH8Pb2hqenJ/vg0Ol0OHHiBN5//32kpKRgwYIFon6IWSIF27hxI/NuWLVqFW7duoUzZ85g9OjROHLkCFQqFW7cuDGUDWrExnaD7pYtWyCXy1mga25uxu7du/XOISKsWLECMpkM+/bt03tOo9Gw9IRCoUBERASysrJMBnAhCQkJSExMRExMjNkfACMJPhgJ3daqq6vh5OSElStXYtasWZg2bZroHYWFq2hTlWvC/DDf0LI/+WGgOxddXl6OCRMmmFzZW8rDhw9RWloKe3t7yOVyVlmXk5ODa9euobGxEW5ubvjwww/Zxp5YWCIFM2fsCMR2g25TUxOWLFnCVlQZGRmQyWRQq9VYvXo1srOzcfnyZcyaNQtBQUHsPwQvDfvZz37GZFNlZWWorq5mQdhUqgLozg/Pnj0bJSUlcHZ2ZqWvPH2NH2l0dnZi9uzZeOGFFxAZGYmvvvoK+fn5+M9//gNHR0dMnz6d5YeVSuWAg9eDBw+Y2c5AVtHCogjeStJY1whhMUVgYKBF7ZlMIbR59Pf310sX6HQ6HD9+HB9//DGmT5+O9vZ2XLt2Ddu3b8fzzz8v2j1YIgWrqqoaDq5gYmO7kjG5XI7z588bPP7UU08hOzsbABAZGcl2TPn8cHJyMpKTk9n5O3fuxIoVK8z2DG1ra8PixYuxb98+A39YPrfb0tKCmJiYXnO7P//5z3Hq1CkoFAqUlJSwx4erdnjUqFE4f/482yibOXMmXn31VRARWlpamAn8jh07UFlZqdetWaVSmfRZALoD0DfffIOGhgYEBASYrJAzhYODA2QymV5w4wsj+B517e3t0Gq1kMlk8PHxsZruljfcGT16NFQqlV4OWq1WIykpCe7u7sjOzh7wfM3BEinYcHcFE5shH3T7y+eff97rcxMmTGDpBo1G06vxh1arxeLFi7Fs2TK2ISccv3fvXoSFhaG6uhrR0dFGtcMA8PLLLyMxMRHLly83eG64aoeNKRM4joOLiwuio6MRHR0N4LvVXX5+PuvWzJcb8/nh6dOnM6mVWq1GTU0NM/UR+yu+k5MTFAoFXF1d8fXXX8POzg7e3t7o6Ohg1XlitRri58/riHsa7uh0Onz88cf4/e9/j1//+teIjY21ujLBEinYcHcFE5thF3RNER8fj2PHjmHr1q04duwYEhISDM4hIqxatQqBgYEGG3L8+KysLCxevBgJCQm9aocBYPbs2aiqqrLWdGwajuPg6ekJT09PLF68GEB3XrG0tBT5+fn4y1/+gm3btuHRo0dwcHDA9773PaSmplotpyo0ZO+ZI+7pYMb3UBO2GuqPHpc3MR83bpyBjvjWrVtISkqCl5cXLl68aDXrzJ5YIgXr7Owc1q5gYjPkc7piYml+mB+fm5uLqKgoNt5UbreqqgoLFy40SC8cPXoUzs7OknbYBPn5+Vi9ejUWLFgAuVyOgoIClJeXY/z48UwpERYWptcVYSBY0qdMqMe9d+8e0+MK9cPCrsK8vnfy5Ml66QI+d3vo0CG88847iImJGdTVoiVSsGHiCiY2truR9rgQSztsLOj2VzsM9J4fHo7aYZ6mpiZ0dXXppYGICLdv39ZTS6jVakycOFHPBN7Z2bnPoGWtPmU9HcE6OjpY+bOzszMCAgL0PCBqamqwceNG+Pr6Yvfu3f3uMdcfTL1fhM5ekZGRyMvLQ3NzMyIjI3Hy5EkQEcLDw1FcXAxHR0d4eXlh3759iIuLGw6uYGIjBV0x6Y90zFjQFWKOdhgALl68iLFjx2L58uV61xrp2mHgOxN43m3typUrrPpNaAIv3Azj2w2NHTsWSqVS1FJhIUSEqqoq1NXVQaFQsJVxfX09jh07BldXVxQUFOC9997Dc889Z/XVrTnvF1MSsJ7etxK9IgVdMTFHO8xjLOgOVDts7FqSdtg4HR0dzAS+sLCQmcAHBwfj7t27cHR0RGpqqlW/FfDVa+PHj4evr69eLrqyshKvv/46OI6DQqFAaWkpwsPD8bvf/c5q9wOY934xJR+Tgq7Z2K5kbCiydetWLFmyBGlpaSw3DEAvNwwAP/3pT5Gbm4vGxkZ4enriV7/6FVatWoXk5GSmHb5z5w4OHjwIwDxviZ7015tipODk5IQZM2YwDw4iwqVLl7B27Vp4eXnByckJsbGx8PDwYKvh0NBQ5k9rCby07fbt2wbVazqdDmlpaThy5Aj27NmDuXPnstfr6uqy6HXNwZz3S18SMMnLxDKsU1ozzOG1wxUVFTh//jzTewq1wwDwySefQKPRYM6cOXB1dcXevXsxdepUFBcXg4jw1ltv6fX3EgOdToeYmBj4+fkhJiam11zzmTNnEBAQAKVSidTUVPb4jh074OHhgZCQEISEhOjNx5bhO/H+7W9/w9mzZ/Hpp5/i2rVrOHToEIKCgnDp0iUsW7YMERERePnll7F//37k5eWhvb29X6/T2tqKoqIiEBFUKpVewK2qqkJ8fDwqKipw+fJlREdH6wV4sVIc8+bNY12phT9ZWVlmjTclAVu/fj0qKytx9epVuLu7Y/PmzaLc80hCWukOAmJoh80db29v36f/cFdXFzZs2KCXs4uPjx/2/sOzZs3SO7azs4OPjw98fHywdOlSAN0abd4E/o9//CO2bNkCOzs7Vk2nUqng7+9vECCF3gw9W+d0dXUhLS0NR48exb59+zBnzhyr5m4tfb+Zko/xnZgBYM2aNay5qoT5SCvdxwyv/QXQq3a4P+OBbt9h/ndmZqbBmIKCAiiVSvj6+mLUqFFYunSp2aug4Y6joyNCQkKwbt06pKWloaioCP/4xz/w4osv4s6dO0hJScGzzz6LhQsXYvv27fj0009x+vRpvPHGG7Czs4NKpdILuDdv3kRCQgJu3ryJL774AlFRUY+1cMCc91tYWBgqKipw8+ZNdHZ24sSJE4iPjwfQvR/Bc/LkSUydOnVwbnw4QUSmfiSsTGNjI82dO5eUSiXNnTuXmpqaiIiotraWnnvuOXbe0qVLyc3NjRwcHMjDw4MOHz5sdLyzs7Pe9V1dXQ1eMyMjg1atWsWOjx8/Ths2bCAiou3bt5O3tzcFBQXRypUrqbm5WfQ52zo6nY5qa2spPT2dIiIiaMKECRQWFkaLFi2it956i3JyckitVtOePXsoJCSELly4QDqdzmr309TURPPmzSOlUknz5s3r9d9s5cqVJJfLacyYMXrvt6amJpo1axaNHj2ajT99+jT5+fmRr68v7dy5k13jpZdeoqlTp1JQUBA9//zzpFarrTYvG6fXuCqpF2wQSzXEGRkZ+Oyzz3D48GEAwB/+8AcUFBRg//79ehriX/7yl9BoNPjoo4+sOyEbJSMjAzU1Ncxvuby8nOmHz5w5g/DwcBw9etRqRjk85soGJdnhoGLand7Ej4SN4e/vz1YfarWa/P39Dc758ssvaf78+ew4JSWFUlJS2HFOTg75+/uTt7c3KRQKg/E6nY42btxIkyZNoqCgILpy5YrB2EmTJtGuXbvEnJpNYc2VbU/M+TfnuXnzJj399NMDHi9hNr3GVSmnO8ywNGd369YtbNiwATk5OdiwYQM6OjpQWlqqNz4nJwcVFRWoqKjABx98gPXr1wP4boMuJycHpaWl+OSTTwzGjhQGM29rqWxQkh0OLpJ6YZhhjobYwcEBBw4cQGxsLCvb5OvkV69ejYaGBixatAg+Pj549dVXkZWVpWdInZWVheXLl4PjOMycORN3796FRqNhvbx8fX0BgG3QjXAza1EwlVKSsC2koDvMMMd/GADi4uKMeviuXr0anp6eevleyVf18TOYssP+jpfoH1J6QUIPEtlXtbq62mgRhvBaSUlJUCqVmDZtGr766iv2nI+PD4KCghASEgKVSjWQ6YwIxJYd9ne8RP+Qgq6E9yo4qQAABA5JREFUHpb4qvZ8vLq6Gl988YXJHG9v+WGef/7zn7h69SqKiorEnKZN0NzcbFZ1YW1tLd544w04OTnh3LlzrKHk5s2b8cQTT7Dqwjlz5iAiIgLl5eXw9PREWloagO6U1Llz5+Dn56c3XsJKmNplG/T9PonHjlarpYkTJ9KNGzeoo6ODpk2bRiUlJXrnnDp1in74wx+STqejf/3rXxQWFmZ07KRJk+jZZ59l43qqJIiIXnnlFfrTn/7EjoU76d7e3nT79m1rTXXIs2XLFqYA2bVrFyUnJxs978KFC3TlyhUDVcL27dvpnXfesfp9ShhFUi9ImIdwky0wMBBLlizB008/jYMHDzJjnri4OPj6+kKpVGLNmjXMGavn2LCwMAQGBrJr87lfIb3lh4HulMX8+fMRGhqKDz74wNpTH3JkZWX1WV0IdHcoEfZ7kxjimIrIj+XzQWLYkJ6eblD5lpiYqHdOXFwcXbp0iR3PnTuXioqKiKi7Ko+IqL6+nqZNm0YXLlzoUwdcVlZGM2fOpFGjRhms8mxNQ+zi4qJ3bKy6kMeY/laqLnysSCtdicHHkvwwAPZboVDgRz/6EfLy8vrUActkMrz//vsGhj1DVUNsqSOYKSRHsKGJFHQlrIapIgye+Ph4HD9+HESEvLw8uLi4wN3dHffv30drayuA7i4PZ8+ehaOjY59GPQqFAmFhYQZ9zoaqyc/nn3+OkpISg5+EhAQm5QIwYCmYvb097OzssGbNGhQUFFhjChL9RAq6ElbDkvxwfX09IiMjERwcjPDwcCxYsABeXl695n/7wlTueKhiqZRLcgQbmkjFERJWxVgRxrp169jfHMfht7/9rcE4X19f/Pvf/9Z7jK+uE2JuuS2ZoT8eaojZocTHxweHDh16nNOR+Ja+XMYkJIYMHMdFANhBRLHfHv8PABDRLiPn7gDQRkTvmhoLoBjAewDsARwmotQe15kM4AiAZwBs46/37XNVAFoBdAF4RERSBYdEn0grXQlbohCAH8dxEwHUAlgK4EULxr4EIBNADIBbAAo5jvs7EQl32JoBJAFY1Mt1f0BEjf2eicSIRQq6EjYDET3iOC4RwGfoXpl+RETXOY5b9+3zBzmOcwNQBMAZgI7juP8GMIWIWnqOBTAWwP8R0Q0A4DjuBIAEAKWC12wA0MBx3IJBm6jEsEYKuhI2BRFlA8ju8dhBwd91ADzNGctx3AsAagSn3AIwoz+3A+Asx3EE4BARjbwKDol+IwVdiZGMsZ20/mxy/BcRqTmOUwA4x3Hc/xLRRZHuTWKYIknGJEYytwB4CY49AajNHUxE6m9/NwA4CSBc1LuTGJZIQVdiJMM21ziOG4XuzbW/mzOQ47gxHMeN4/8GMB9AielREhLA/wMq4qQ84McfWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "x1 = all_data.iloc[:,1]\n",
    "x2 = all_data.iloc[:,2]\n",
    "x3 = all_data.iloc[:,3]\n",
    "ax.scatter3D(x1, x2, x3, c=x3, cmap='Greens');\n",
    "\n",
    "plt.show()\n",
    "\n",
    "#actual implementation of liner regression\n",
    "#compute y_pred, compare with y_data above etc etc\n",
    "#assume more code here\n",
    "\n",
    "exit() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept     0.001327\n",
      "mkt_excess    0.897902\n",
      "SMB           0.027838\n",
      "HML          -0.113354\n",
      "dtype: float64\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:            port_excess   R-squared:                       0.798\n",
      "Model:                            OLS   Adj. R-squared:                  0.797\n",
      "Method:                 Least Squares   F-statistic:                     636.7\n",
      "Date:                Wed, 18 Nov 2020   Prob (F-statistic):          1.50e-167\n",
      "Time:                        15:12:58   Log-Likelihood:                 1196.6\n",
      "No. Observations:                 488   AIC:                            -2385.\n",
      "Df Residuals:                     484   BIC:                            -2368.\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept      0.0013      0.001      1.380      0.168      -0.001       0.003\n",
      "mkt_excess     0.8979      0.022     40.577      0.000       0.854       0.941\n",
      "SMB            0.0278      0.033      0.831      0.406      -0.038       0.094\n",
      "HML           -0.1134      0.033     -3.445      0.001      -0.178      -0.049\n",
      "==============================================================================\n",
      "Omnibus:                      238.735   Durbin-Watson:                   2.152\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             2663.014\n",
      "Skew:                          -1.838   Prob(JB):                         0.00\n",
      "Kurtosis:                      13.837   Cond. No.                         38.3\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "0.1683525515061144\n",
      "7.151419559919375e-158\n",
      "0.40636552386776614\n",
      "0.0006200078819804644\n",
      "1.379581549021653\n",
      "40.577247797576206\n",
      "0.8310316375383785\n",
      "-3.445320444601942\n",
      "0.0013265239002386762\n"
     ]
    }
   ],
   "source": [
    "model = sm.formula.ols(formula = \"port_excess ~ mkt_excess + SMB + HML\", data = all_data).fit()\n",
    "print(model.params)\n",
    "print(model.summary())\n",
    "print(model.pvalues[0]) # alpha\n",
    "print(model.pvalues[1]) # market excess\n",
    "print(model.pvalues[2]) # SMB\n",
    "print(model.pvalues[3]) # HML\n",
    "#model??\n",
    "#display(help(model))\n",
    "print(model.tvalues[0]) # alpha\n",
    "print(model.tvalues[1]) # market excess\n",
    "print(model.tvalues[2]) # SMB\n",
    "print(model.tvalues[3]) # HML\n",
    "print(model.params.Intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_reg_model(ticker,start,end):\n",
    "    # Get FF data\n",
    "    ff_data = get_fama_french()\n",
    "    ff_last = ff_data.index[ff_data.shape[0] - 1].date()\n",
    "    #Get the fund price data\n",
    "    price_data = get_price_data(ticker,start,end)\n",
    "    price_data = price_data.loc[:ff_last]\n",
    "    ret_data = get_return_data(price_data, \"M\")\n",
    "    all_data = pd.merge(pd.DataFrame(ret_data),ff_data, how = 'inner', left_index= True, right_index= True)\n",
    "    all_data.rename(columns={\"Mkt-RF\":\"mkt_excess\"}, inplace=True)\n",
    "    all_data['port_excess'] = all_data['portfolio'] - all_data['RF']\n",
    "    # Run the model\n",
    "    model = sm.formula.ols(formula = \"port_excess ~ mkt_excess + SMB + HML\", data = all_data).fit()\n",
    "    #return model.params\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<statsmodels.regression.linear_model.RegressionResultsWrapper object at 0x000001BA48038E80>\n",
      "1.472300617687898e-12\n"
     ]
    }
   ],
   "source": [
    "ggrax_model = run_reg_model(\"GGRAX\", start = \"1999-05-01\", end = \"2019-06-30\")\n",
    "print(ggrax_model)\n",
    "print(ggrax_model.pvalues[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3849\n",
      "1888\n",
      "831\n",
      "234\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Symbol</th>\n",
       "      <th>Name</th>\n",
       "      <th>LastSale</th>\n",
       "      <th>MarketCap</th>\n",
       "      <th>IPOyear</th>\n",
       "      <th>Sector</th>\n",
       "      <th>industry</th>\n",
       "      <th>Summary Quote</th>\n",
       "      <th>Unnamed: 8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>FLWS</td>\n",
       "      <td>1-800-FLOWERS.COM, Inc.</td>\n",
       "      <td>21.37</td>\n",
       "      <td>$1.38B</td>\n",
       "      <td>1999.0</td>\n",
       "      <td>Consumer Services</td>\n",
       "      <td>Other Specialty Stores</td>\n",
       "      <td>https://old.nasdaq.com/symbol/flws</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>JOBS</td>\n",
       "      <td>51job, Inc.</td>\n",
       "      <td>69.70</td>\n",
       "      <td>$4.68B</td>\n",
       "      <td>2004.0</td>\n",
       "      <td>Technology</td>\n",
       "      <td>Diversified Commercial Services</td>\n",
       "      <td>https://old.nasdaq.com/symbol/jobs</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>ACAD</td>\n",
       "      <td>ACADIA Pharmaceuticals Inc.</td>\n",
       "      <td>47.35</td>\n",
       "      <td>$7.52B</td>\n",
       "      <td>1985.0</td>\n",
       "      <td>Health Care</td>\n",
       "      <td>Major Pharmaceuticals</td>\n",
       "      <td>https://old.nasdaq.com/symbol/acad</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>72</th>\n",
       "      <td>ADBE</td>\n",
       "      <td>Adobe Inc.</td>\n",
       "      <td>461.11</td>\n",
       "      <td>$221.2B</td>\n",
       "      <td>1986.0</td>\n",
       "      <td>Technology</td>\n",
       "      <td>Computer Software: Prepackaged Software</td>\n",
       "      <td>https://old.nasdaq.com/symbol/adbe</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75</th>\n",
       "      <td>AEIS</td>\n",
       "      <td>Advanced Energy Industries, Inc.</td>\n",
       "      <td>69.02</td>\n",
       "      <td>$2.65B</td>\n",
       "      <td>1995.0</td>\n",
       "      <td>Capital Goods</td>\n",
       "      <td>Industrial Machinery/Components</td>\n",
       "      <td>https://old.nasdaq.com/symbol/aeis</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Symbol                              Name  LastSale MarketCap  IPOyear  \\\n",
       "7    FLWS           1-800-FLOWERS.COM, Inc.     21.37    $1.38B   1999.0   \n",
       "17   JOBS                       51job, Inc.     69.70    $4.68B   2004.0   \n",
       "33   ACAD       ACADIA Pharmaceuticals Inc.     47.35    $7.52B   1985.0   \n",
       "72   ADBE                        Adobe Inc.    461.11   $221.2B   1986.0   \n",
       "75   AEIS  Advanced Energy Industries, Inc.     69.02    $2.65B   1995.0   \n",
       "\n",
       "               Sector                                 industry  \\\n",
       "7   Consumer Services                   Other Specialty Stores   \n",
       "17         Technology          Diversified Commercial Services   \n",
       "33        Health Care                    Major Pharmaceuticals   \n",
       "72         Technology  Computer Software: Prepackaged Software   \n",
       "75      Capital Goods          Industrial Machinery/Components   \n",
       "\n",
       "                         Summary Quote  Unnamed: 8  \n",
       "7   https://old.nasdaq.com/symbol/flws         NaN  \n",
       "17  https://old.nasdaq.com/symbol/jobs         NaN  \n",
       "33  https://old.nasdaq.com/symbol/acad         NaN  \n",
       "72  https://old.nasdaq.com/symbol/adbe         NaN  \n",
       "75  https://old.nasdaq.com/symbol/aeis         NaN  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 41.8 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['FLWS',\n",
       " 'JOBS',\n",
       " 'ACAD',\n",
       " 'ADBE',\n",
       " 'AEIS',\n",
       " 'AKAM',\n",
       " 'ALXN',\n",
       " 'ALGN',\n",
       " 'ALKS',\n",
       " 'ALGT',\n",
       " 'ALNY',\n",
       " 'GOOG',\n",
       " 'AIMC',\n",
       " 'AMZN',\n",
       " 'AMWD',\n",
       " 'ABCB',\n",
       " 'AMSF',\n",
       " 'AMGN',\n",
       " 'AMKR',\n",
       " 'ANSS',\n",
       " 'AAPL',\n",
       " 'AMAT',\n",
       " 'ARNA',\n",
       " 'ASML',\n",
       " 'AZPN',\n",
       " 'ASTE',\n",
       " 'ATRC',\n",
       " 'BIDU',\n",
       " 'BANF',\n",
       " 'BECN',\n",
       " 'BBBY',\n",
       " 'BMRN',\n",
       " 'BLKB',\n",
       " 'EPAY',\n",
       " 'BRKS',\n",
       " 'BLDR',\n",
       " 'CHRW',\n",
       " 'CALM',\n",
       " 'CSIQ',\n",
       " 'CSII',\n",
       " 'CWST',\n",
       " 'CASY',\n",
       " 'CENT',\n",
       " 'CERN',\n",
       " 'GTLS',\n",
       " 'CHKP',\n",
       " 'CTAS',\n",
       " 'CRUS',\n",
       " 'CSCO',\n",
       " 'CTXS',\n",
       " 'CCMP',\n",
       " 'CME',\n",
       " 'COKE',\n",
       " 'CGNX',\n",
       " 'CTSH',\n",
       " 'COLB',\n",
       " 'COLM',\n",
       " 'CVLT',\n",
       " 'CGEN',\n",
       " 'CNMD',\n",
       " 'CPRT',\n",
       " 'CORT',\n",
       " 'CSGP',\n",
       " 'CBRL',\n",
       " 'CACC',\n",
       " 'CREE',\n",
       " 'CROX',\n",
       " 'CSGS',\n",
       " 'CYTK',\n",
       " 'XRAY',\n",
       " 'DXCM',\n",
       " 'DISH',\n",
       " 'DLTR',\n",
       " 'EBAY',\n",
       " 'EHTH',\n",
       " 'ENTG',\n",
       " 'EQIX',\n",
       " 'EEFT',\n",
       " 'EXAS',\n",
       " 'EXEL',\n",
       " 'EXLS',\n",
       " 'FFIV',\n",
       " 'FARO',\n",
       " 'FAST',\n",
       " 'FSLR',\n",
       " 'FCFS',\n",
       " 'FISV',\n",
       " 'FLEX',\n",
       " 'FLIR',\n",
       " 'FORM',\n",
       " 'GRMN',\n",
       " 'ROCK',\n",
       " 'GILD',\n",
       " 'OMAB',\n",
       " 'HCSG',\n",
       " 'HTLD',\n",
       " 'HSIC',\n",
       " 'HSKA',\n",
       " 'HMSY',\n",
       " 'HOLX',\n",
       " 'HOMB',\n",
       " 'HUBG',\n",
       " 'ICFI',\n",
       " 'ICLR',\n",
       " 'ICUI',\n",
       " 'IDXX',\n",
       " 'INFO',\n",
       " 'IIVI',\n",
       " 'ILMN',\n",
       " 'IMGN',\n",
       " 'NSIT',\n",
       " 'INTU',\n",
       " 'ISRG',\n",
       " 'IPGP',\n",
       " 'IRBT',\n",
       " 'ITRI',\n",
       " 'JJSF',\n",
       " 'JCOM',\n",
       " 'JKHY',\n",
       " 'JBLU',\n",
       " 'KLAC',\n",
       " 'LRCX',\n",
       " 'LAMR',\n",
       " 'LSTR',\n",
       " 'LHCG',\n",
       " 'LPSN',\n",
       " 'LOGI',\n",
       " 'LMNX',\n",
       " 'MANH',\n",
       " 'MANT',\n",
       " 'MKTX',\n",
       " 'MRTN',\n",
       " 'MRVL',\n",
       " 'MXIM',\n",
       " 'MGRC',\n",
       " 'MRCY',\n",
       " 'CASH',\n",
       " 'MCHP',\n",
       " 'MSFT',\n",
       " 'MSTR',\n",
       " 'MKSI',\n",
       " 'MPWR',\n",
       " 'MNRO',\n",
       " 'MORN',\n",
       " 'NATI',\n",
       " 'NEOG',\n",
       " 'NTAP',\n",
       " 'NFLX',\n",
       " 'NTCT',\n",
       " 'NBIX',\n",
       " 'NXST',\n",
       " 'EGOV',\n",
       " 'NVMI',\n",
       " 'NUVA',\n",
       " 'NVDA',\n",
       " 'ODFL',\n",
       " 'OMCL',\n",
       " 'OTEX',\n",
       " 'ORLY',\n",
       " 'OSIS',\n",
       " 'OSTK',\n",
       " 'PZZA',\n",
       " 'PDCO',\n",
       " 'PAYX',\n",
       " 'PEGA',\n",
       " 'PENN',\n",
       " 'PRFT',\n",
       " 'PLUG',\n",
       " 'POOL',\n",
       " 'POWI',\n",
       " 'PRAA',\n",
       " 'PFG',\n",
       " 'PRGS',\n",
       " 'QCOM',\n",
       " 'RDWR',\n",
       " 'RMBS',\n",
       " 'ROLL',\n",
       " 'REG',\n",
       " 'REGN',\n",
       " 'RCII',\n",
       " 'RGEN',\n",
       " 'ROST',\n",
       " 'RYAAY',\n",
       " 'SAFT',\n",
       " 'SAFM',\n",
       " 'SGMO',\n",
       " 'SANM',\n",
       " 'SPNS',\n",
       " 'SBAC',\n",
       " 'STX',\n",
       " 'SGEN',\n",
       " 'SEIC',\n",
       " 'SBNY',\n",
       " 'SLGN',\n",
       " 'SLAB',\n",
       " 'SIMO',\n",
       " 'SBGI',\n",
       " 'SKYW',\n",
       " 'STMP',\n",
       " 'SBUX',\n",
       " 'STLD',\n",
       " 'SRCL',\n",
       " 'SHOO',\n",
       " 'STRA',\n",
       " 'SPWR',\n",
       " 'SYKE',\n",
       " 'SYNA',\n",
       " 'SNPS',\n",
       " 'TROW',\n",
       " 'TTWO',\n",
       " 'TTEK',\n",
       " 'TCBI',\n",
       " 'TXRH',\n",
       " 'CAKE',\n",
       " 'PRSC',\n",
       " 'TSEM',\n",
       " 'TSCO',\n",
       " 'TBIO',\n",
       " 'TRMB',\n",
       " 'TTEC',\n",
       " 'TTMI',\n",
       " 'UFPI',\n",
       " 'UTHR',\n",
       " 'URBN',\n",
       " 'VRNT',\n",
       " 'VRSN',\n",
       " 'VRTX',\n",
       " 'VSAT',\n",
       " 'VG',\n",
       " 'WDFC',\n",
       " 'WERN',\n",
       " 'WYNN',\n",
       " 'XLNX',\n",
       " 'ZBRA']"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "tickers_df = pd.read_csv('companylist.csv')\n",
    "print(len(tickers_df))\n",
    "tickers_df = tickers_df.loc[tickers_df['MarketCap'].str.contains(r'B',na=True)]\n",
    "print(len(tickers_df))\n",
    "tickers_df = tickers_df.dropna(subset=['Sector','industry','LastSale','IPOyear'],axis=0,how='any')\n",
    "print(len(tickers_df))\n",
    "tickers_df = tickers_df.loc[tickers_df['IPOyear']<2007]\n",
    "print(len(tickers_df))\n",
    "display(tickers_df.head())\n",
    "tickers_list = tickers_df.iloc[:,0].tolist()\n",
    "tickers_list\n",
    "# tickers_str = \" \".join(tickers_list)\n",
    "# tickers_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "233\n",
      "232\n",
      "231\n",
      "230\n",
      "229\n",
      "228\n",
      "227\n",
      "226\n",
      "225\n",
      "224\n",
      "223\n",
      "222\n",
      "221\n",
      "220\n",
      "219\n",
      "218\n",
      "217\n",
      "216\n",
      "215\n",
      "214\n",
      "213\n",
      "212\n",
      "211\n",
      "210\n",
      "209\n",
      "208\n",
      "207\n",
      "206\n",
      "205\n",
      "204\n",
      "203\n",
      "202\n",
      "201\n",
      "200\n",
      "199\n",
      "198\n",
      "197\n",
      "196\n",
      "195\n",
      "194\n",
      "193\n",
      "192\n",
      "191\n",
      "190\n",
      "189\n",
      "188\n",
      "187\n",
      "186\n",
      "185\n",
      "184\n",
      "183\n",
      "182\n",
      "181\n",
      "180\n",
      "179\n",
      "178\n",
      "177\n",
      "176\n",
      "175\n",
      "174\n",
      "173\n",
      "172\n",
      "171\n",
      "170\n",
      "169\n",
      "168\n",
      "167\n",
      "166\n",
      "165\n",
      "164\n",
      "163\n",
      "162\n",
      "161\n",
      "160\n",
      "159\n",
      "158\n",
      "157\n",
      "156\n",
      "155\n",
      "154\n",
      "153\n",
      "152\n",
      "151\n",
      "150\n",
      "149\n",
      "148\n",
      "147\n",
      "146\n",
      "145\n",
      "144\n",
      "143\n",
      "142\n",
      "141\n",
      "140\n",
      "139\n",
      "138\n",
      "137\n",
      "136\n",
      "135\n",
      "134\n",
      "133\n",
      "132\n",
      "131\n",
      "130\n",
      "129\n",
      "128\n",
      "127\n",
      "126\n",
      "125\n",
      "124\n",
      "123\n",
      "122\n",
      "121\n",
      "120\n",
      "119\n",
      "118\n",
      "117\n",
      "116\n",
      "115\n",
      "114\n",
      "113\n",
      "112\n",
      "111\n",
      "110\n",
      "109\n",
      "108\n",
      "107\n",
      "106\n",
      "105\n",
      "104\n",
      "103\n",
      "102\n",
      "101\n",
      "100\n",
      "99\n",
      "98\n",
      "97\n",
      "96\n",
      "95\n",
      "94\n",
      "93\n",
      "92\n",
      "91\n",
      "90\n",
      "89\n",
      "88\n",
      "87\n",
      "86\n",
      "85\n",
      "84\n",
      "83\n",
      "82\n",
      "81\n",
      "80\n",
      "79\n",
      "78\n",
      "77\n",
      "76\n",
      "75\n",
      "74\n",
      "73\n",
      "72\n",
      "71\n",
      "70\n",
      "69\n",
      "68\n",
      "67\n",
      "66\n",
      "65\n",
      "64\n",
      "63\n",
      "62\n",
      "61\n",
      "60\n",
      "59\n",
      "58\n",
      "57\n",
      "56\n",
      "55\n",
      "54\n",
      "53\n",
      "52\n",
      "51\n",
      "50\n",
      "49\n",
      "48\n",
      "47\n",
      "46\n",
      "45\n",
      "44\n",
      "43\n",
      "42\n",
      "41\n",
      "40\n",
      "39\n",
      "38\n",
      "37\n",
      "36\n",
      "35\n",
      "34\n",
      "33\n",
      "32\n",
      "31\n",
      "30\n",
      "29\n",
      "28\n",
      "27\n",
      "26\n",
      "25\n",
      "24\n",
      "23\n",
      "22\n",
      "21\n",
      "20\n",
      "19\n",
      "18\n",
      "17\n",
      "16\n",
      "15\n",
      "14\n",
      "13\n",
      "12\n",
      "11\n",
      "10\n",
      "9\n",
      "8\n",
      "7\n",
      "6\n",
      "5\n",
      "4\n",
      "3\n",
      "2\n",
      "1\n",
      "0\n",
      "0.04700854700854701%\n",
      "Wall time: 14min 40s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['ALXN',\n",
       " 'CHKP',\n",
       " 'GILD',\n",
       " 'HCSG',\n",
       " 'HOMB',\n",
       " 'ILMN',\n",
       " 'MNRO',\n",
       " 'NBIX',\n",
       " 'NVDA',\n",
       " 'ODFL',\n",
       " 'POOL']"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "a=0\n",
    "b=0\n",
    "ticker=[]\n",
    "for i in tickers_list[:]:\n",
    "    try:\n",
    "        ggrax_model = run_reg_model(i, start = \"1999-05-01\", end = \"2019-06-30\")\n",
    "        a+=1\n",
    "        if ggrax_model.pvalues[0]<0.05 and ggrax_model.pvalues[2]<0.05 and ggrax_model.pvalues[3]<0.05 and ggrax_model.params.Intercept>0:\n",
    "            b+=1\n",
    "            ticker.append(i)\n",
    "        print(234-a) #234 after filtering nasdaq stocks\n",
    "    except Exception as e:\n",
    "        print(str(e))\n",
    "        pass\n",
    "print(str(b/a)+'%')\n",
    "ticker #statistically significant p value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.019564266466639814\n",
      "0.013336366276620437\n",
      "0.01595993235517934\n",
      "0.011226784488733519\n",
      "0.010006630750737407\n",
      "0.020318938219351008\n",
      "0.0133587832870867\n",
      "0.022563840944168086\n",
      "0.023828074552818386\n",
      "0.016747473392872927\n",
      "0.012958578778308109\n"
     ]
    }
   ],
   "source": [
    "for i in ['ALXN','CHKP','GILD','HCSG','HOMB','ILMN','MNRO','NBIX','NVDA','ODFL','POOL']:\n",
    "    ggrax_model = run_reg_model(i, start = \"1999-05-01\", end = \"2019-06-30\")\n",
    "    print(ggrax_model.params.Intercept)\n",
    "\n",
    "#  find highest 3 R^2 stocks since we only have 3 unknowns:\n",
    "# .335 ALXN -65.87% alpha=0.019566\n",
    "# .361 CHKP 111.3% alpha=0.01333\n",
    "# .174\n",
    "# .121\n",
    "# .378 HOMB 54.57% alpha=0.01\n",
    "# .271\n",
    "# .151\n",
    "# .260\n",
    "# .314\n",
    "# .155\n",
    "# .319\n",
    "          \n",
    "ticker = ['ALXN','CHKP','HOMB']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.92758213  1.07307085  0.56230307]\n",
      " [ 1.61951529  0.60609821  0.72037254]\n",
      " [-1.50440748 -1.35598889  0.94509794]]\n",
      "[[-1 -1 -1]\n",
      " [ 0  0  0]\n",
      " [ 0  0  0]]\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "beta=np.array([])\n",
    "smb=np.array([])\n",
    "hml=np.array([])\n",
    "for i in ticker[:3]:\n",
    "    beta = np.append(beta,[run_reg_model(i, start = \"1999-05-01\", end = \"2019-06-30\").params.mkt_excess])\n",
    "    smb = np.append(smb,[run_reg_model(i, start = \"1999-05-01\", end = \"2019-06-30\").params.SMB])\n",
    "    hml = np.append(hml,[run_reg_model(i, start = \"1999-05-01\", end = \"2019-06-30\").params.HML])\n",
    "# display(beta)\n",
    "# display(smb)\n",
    "# display(hml)\n",
    "matrix = np.vstack([beta,smb,hml])\n",
    "print(matrix) #number of asset by 3 independent variables\n",
    "#z = np.zeros([3,3])\n",
    "z = np.array([[-1]*3,[0]*3,[0]*3])\n",
    "print(z)\n",
    "\n",
    "#matrix = np.append(matrix,z,axis=1) #add column of zeros to form square matrix\n",
    "#matrix = np.linalg.solve(matrix,z) # only for square matrix\n",
    "np.linalg.inv(matrix).dot(z) # only for square matrix\n",
    "#https://scipy.github.io/devdocs/tutorial/linalg.html\n",
    "\n",
    "#x, residuals, rank, s = np.linalg.lstsq(matrix,z) #Least squares is a standard approach to problems with more equations than unknowns, also known as overdetermined systems.\n",
    "#x, residuals, rank, singular version of matrix(also absolute of eigenvalue)\n",
    "#https://sodocumentation.net/numpy/topic/3753/linear-algebra-with-np-linalg\n",
    "print(x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # optimization module\n",
    "# def objective(x): # x are the p-value significant nasdaq stocks\n",
    "#     x1=x[0]\n",
    "#     x2=x[1]\n",
    "#     x3=x[2]\n",
    "#     x4=x[3]\n",
    "#     return -1.50440748*x1-1.35598889*x2-0.50999618*x3+0.77768992*x4\n",
    "# def constraint1(x):\n",
    "#     return -x[0]-x[1]-x[2]-x[3]\n",
    "# initializer = np.ones([4,1])\n",
    "# #print(objective(initializer))\n",
    "# bounds = tuple((-1,1) for i in range(0,4))\n",
    "# con1 = {'type':'eq','fun':constraint1}\n",
    "# cons = [con1]\n",
    "# #print(initializer)\n",
    "# #print(bounds)\n",
    "\n",
    "# sol = minimize(objective,initializer,method='SLSQP',bounds=bounds,constraints=[])\n",
    "# display(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}