Projet_MR_2023_final.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projet Market Risk 2023\n",
    "\n",
    "- Clément Muret\n",
    "- Maxime Le Floch\n",
    "- IF3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_path = r\"C:\\Users\\clemm\\OneDrive\\Documents\\Cours\\S7\\Market risk\\Natixis stock.txt\"\n",
    "\n",
    "# Dictionary for storing prices by date\n",
    "price_data = {}\n",
    "\n",
    "with open(file_path, 'r') as file:\n",
    "    lines = file.readlines()\n",
    "\n",
    "for line in lines:\n",
    "    date, price = line.strip().split('\\t')  # Tab split\n",
    "    # Conversions\n",
    "    from datetime import datetime\n",
    "    date = datetime.strptime(date, \"%d/%m/%Y\")\n",
    "    price = float(price.replace(',', '.'))\n",
    "\n",
    "    # Storing prices in the dictionary by date\n",
    "    if date.year not in price_data:\n",
    "        price_data[date.year] = {}\n",
    "    if date.month not in price_data[date.year]:\n",
    "        price_data[date.year][date.month] = []\n",
    "\n",
    "    price_data[date.year][date.month].append(price)\n",
    "\n",
    "\n",
    "# Retrieve prices from January 2015 to December 2016\n",
    "Prix_2015_2016 = []\n",
    "for year in range(2015, 2017):\n",
    "    for month in range(1, 13):\n",
    "        Prix_2015_2016 += price_data[year][month]\n",
    "\n",
    "# Retrieve prices from January 2017 to December 2018\n",
    "Prix_2017_2018 = []\n",
    "for year in range(2017, 2019):\n",
    "    for month in range(1, 13):\n",
    "        Prix_2017_2018 += price_data[year][month]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Useful functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_returns(prices):\n",
    "    returns = []\n",
    "    for i in range(1, len(prices)):\n",
    "        daily_return = (prices[i] - prices[i - 1]) / prices[i - 1]\n",
    "        returns.append(daily_return)\n",
    "    return returns\n",
    "\n",
    "def calculate_mean(data):\n",
    "    return sum(data) / len(data)\n",
    "\n",
    "def calculate_std_dev(data):\n",
    "    mean = calculate_mean(data)\n",
    "    variance = sum((x - mean) ** 2 for x in data) / len(data)\n",
    "    return variance ** 0.5\n",
    "\n",
    "def calculate_quantile(data, p):\n",
    "    sorted_data = sorted(data)\n",
    "    index = int(p * len(sorted_data))\n",
    "    return sorted_data[index]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question A (Ex2, part of Q1 and of Q2 of TD1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated VaR for 2015-2016 :\n",
      "Confiance\tVaR_loss\t\tVaR_gain\n",
      "90.0%\t\t-0.02780410742496053\t0.02819865319865327\n",
      "95.0%\t\t-0.038933559035124875\t0.03788316946211676\n",
      "99.0%\t\t-0.05704225352112663\t0.055596196049743904\n"
     ]
    }
   ],
   "source": [
    "# This part is here to check that the results found with the kernel and EVT are consistent\n",
    "\n",
    "# Return calculation  \n",
    "returns_2015_2016 = calculate_returns(Prix_2015_2016)\n",
    "\n",
    "confidence_level = [0.1, 0.05, 0.01]\n",
    "print(\"Estimated VaR for 2015-2016 :\")\n",
    "print(\"Confiance\\tVaR_loss\\t\\tVaR_gain\")\n",
    "for p in confidence_level:\n",
    "    loss = calculate_quantile(returns_2015_2016, p)\n",
    "    gain = calculate_quantile(returns_2015_2016, 1 - p)\n",
    "    print(f\"{(1-p)*100}%\\t\\t{loss}\\t{gain}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def biweight_kernel(u):\n",
    "    if abs(u) <= 1:\n",
    "        return (15 / 16) * (1 - u ** 2) ** 2\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "def kernel_density_estimate(x, data, h):\n",
    "    n = len(data)\n",
    "    total = 0\n",
    "    for xi in data:\n",
    "        total += biweight_kernel((x - xi) / h)\n",
    "    return total / (n * h)\n",
    "\n",
    "# Estimate an integral by trapzes method\n",
    "def trapezes(f, a, b, n, data, h):\n",
    "    step = (b - a) / n\n",
    "    sum = 0\n",
    "    for i in range(1, n):\n",
    "        xi = a + i * step\n",
    "        xj = a + (i + 1) * step\n",
    "        sum += f(xi, data, h) + f(xj, data, h)\n",
    "    return sum * step / 2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGzCAYAAACPa3XZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAADMWUlEQVR4nOzdd3hUZdrA4d+Zycyk914IoZPQe0CKdEQQGyoqYllXxUUW21o+xd57XUUFV7CgYlc6SJMSCL2GBEJCejLpkynn+2OSgZgECCQZkjz3dc2VzDnve85zEkievFVRVVVFCCGEEKKJaJwdgBBCCCFaF0k+hBBCCNGkJPkQQgghRJOS5EMIIYQQTUqSDyGEEEI0KUk+hBBCCNGkJPkQQgghRJOS5EMIIYQQTUqSDyGEEEI0KUk+hGhGRowYwYgRI867vqIozJ07t8HiOV9t27ZlxowZzg6jQaSkpKAoCvPnz3d2KEI0G5J8iFZl/vz5KIrCtm3bqh03Go0MGDAAV1dX/vjjDydF13rt27ePuXPnkpKS4uxQ6rRo0SLefPNNZ4dRTXp6OnPnziUxMdHZoQhRLy7ODkAIZyssLGTs2LHs2rWLJUuWMH78eGeH1OIdPHgQjebU3z779u3jqaeeYsSIEbRt29Z5gZ3BokWL2LNnD7Nnz652PDo6mrKyMnQ6XZPHlJ6ezlNPPUXbtm3p1atXk99fiPMlyYdo1YqKihg3bhyJiYl8//33TJgw4YKvWV5ejl6vr/bLVVRnMBicHUKDURQFV1dXZ4chRLMiPx1Fq1VcXMz48ePZvn073333HRMnTqx2Pi0tjdtuu42QkBAMBgNxcXF8+umn1cqsWbMGRVH46quvePzxx4mIiMDd3Z3CwkJmzJiBp6cnaWlpTJkyBU9PT4KCgnjggQewWq3VrmOz2XjzzTeJi4vD1dWVkJAQ/vnPf5Kfn39ez2Yymfj3v/9NUFAQXl5eTJ48mRMnTtRatj7P+c033/Dcc88RGRmJq6sro0aN4siRI9XKHj58mKuvvprQ0FBcXV2JjIzk+uuvx2g0OsqcPuZj/vz5XHvttQBceumlKIqCoiisWbOGW265hcDAQMxmc424x44dS+fOnc/6tdi8eTPjx4/Hx8cHd3d3hg8fzoYNG6qVKSoqYvbs2bRt2xaDwUBwcDBjxoxh+/btgH2sza+//sqxY8cc8VW10NQ25qPqe3/8+HEuv/xyPD09iYiI4L333gNg9+7djBw5Eg8PD6Kjo1m0aFG1ePLy8njggQfo3r07np6eeHt7M2HCBHbu3Fnte9K/f38Abr31Vkdcp8dxLs8uhDNIy4dolUpKSpgwYQJbt27l22+/5fLLL692PjMzk0GDBqEoCvfeey9BQUH8/vvv3H777RQWFtZoen/mmWfQ6/U88MADmEwm9Ho9AFarlXHjxjFw4EBeffVVVqxYwWuvvUb79u25++67HfX/+c9/Mn/+fG699VZmzZpFcnIy7777Ljt27GDDhg31btK/4447+OKLL5g2bRqDBw9m1apVNZKr83nOF198EY1GwwMPPIDRaOTll1/mxhtvZPPmzQBUVFQwbtw4TCYT//rXvwgNDSUtLY1ffvmFgoICfHx8asQwbNgwZs2axdtvv82jjz5K165dAejatSs333wzn3/+OUuXLq32PcrIyGDVqlU8+eSTZ/w6rFq1igkTJtC3b1+efPJJNBoNn332GSNHjmTdunUMGDAAgLvuuotvv/2We++9l9jYWHJzc1m/fj379++nT58+PPbYYxiNRk6cOMEbb7wBgKen5xnvbbVamTBhAsOGDePll19m4cKF3HvvvXh4ePDYY49x4403ctVVV/Hhhx8yffp04uPjiYmJAeDo0aP88MMPXHvttcTExJCZmcl///tfhg8fzr59+wgPD6dr1648/fTTPPHEE9x5550MHToUgMGDB9fr2YVwClWIVuSzzz5TATU6OlrV6XTqDz/8UGu522+/XQ0LC1NzcnKqHb/++utVHx8ftbS0VFVVVV29erUKqO3atXMcq3LLLbeogPr0009XO967d2+1b9++jvfr1q1TAXXhwoXVyv3xxx81jg8fPlwdPnz4GZ8xMTFRBdR77rmn2vFp06apgPrkk0+e93N27dpVNZlMjnJvvfWWCqi7d+9WVVVVd+zYoQLq4sWLzxhjdHS0essttzjeL168WAXU1atXVytntVrVyMhI9brrrqt2/PXXX1cVRVGPHj1a5z1sNpvasWNHddy4carNZnMcLy0tVWNiYtQxY8Y4jvn4+KgzZ848Y8wTJ05Uo6OjaxxPTk5WAfWzzz5zHKv63j///POOY/n5+aqbm5uqKIr61VdfOY4fOHCgxvelvLxctVqtNe5jMBiq/XvaunVrjXvX99mFcAbpdhGtUmZmJq6urkRFRdU4p6oq3333HZMmTUJVVXJychyvcePGYTQaHc3xVW655Rbc3Nxqvdddd91V7f3QoUM5evSo4/3ixYvx8fFhzJgx1e7Vt29fPD09Wb16db2e7bfffgNg1qxZ1Y7/vRXjfJ7z1ltvdbTqVD0L4HieqpaNpUuXUlpaWq+4a6PRaLjxxhv56aefKCoqchxfuHAhgwcPdrQU1CYxMZHDhw8zbdo0cnNzHc9WUlLCqFGj+PPPP7HZbAD4+vqyefNm0tPTLzjm091xxx2Oz319fencuTMeHh5MnTrVcbxz5874+vpW+zdhMBgcY4asViu5ubl4enrSuXPnGt+T2tTn2YVwBkk+RKv03//+F71ez/jx4zl48GC1c9nZ2RQUFPDRRx8RFBRU7XXrrbcCkJWVVa1OXb8EXV1dCQoKqnbMz8+v2liOw4cPYzQaCQ4OrnG/4uLiGvc6m2PHjqHRaGjfvn21438fH3E+z9mmTZsazwI4nicmJoY5c+Ywb948AgMDGTduHO+991618R71NX36dMrKyliyZAlgnymTkJDAzTfffMZ6hw8fBuyJ4d+fb968eZhMJkdcL7/8Mnv27CEqKooBAwYwd+7casnA+ajte+/j40NkZCSKotQ4fvq/CZvNxhtvvEHHjh0xGAwEBgYSFBTErl27zulrWZ9nF8IZZMyHaJViY2P57bffGDVqFGPGjGHDhg2OVpCqvwhvuukmbrnlllrr9+jRo9r7ulo9tFrtWWOx2WwEBwezcOHCWs///RdYQzmf56zreVRVdXz+2muvMWPGDH788UeWLVvGrFmzeOGFF/jrr7+IjIysd5yxsbH07duXL774gunTp/PFF1+g1+urtR7Upur5XnnllTqnoVaN25g6dSpDhw5lyZIlLFu2jFdeeYWXXnrpgmZA1fW1Opev4fPPP8///d//cdttt/HMM8/g7++PRqNh9uzZ59RiUZ9nF8IZJPkQrdaAAQP44YcfmDhxImPGjGHdunWOvw69vLywWq2MHj260eNo3749K1asYMiQIXUmMfURHR2NzWYjKSmpWmvH31t4GvM5u3fvTvfu3Xn88cfZuHEjQ4YM4cMPP+TZZ5+ttfzfWwL+bvr06cyZM4eTJ0+yaNEiJk6c6Gh1qUtVy4+3t/c5PV9YWBj33HMP99xzD1lZWfTp04fnnnvOkXycLcaG9O2333LppZfyySefVDteUFBAYGCg431dMdX32YVoatLtIlq1UaNG8eWXX3LkyBHGjx9PYWEhWq2Wq6++mu+++449e/bUqJOdnd2gMUydOhWr1cozzzxT45zFYqGgoKBe16v6Zfn2229XO/731Tkb4zkLCwuxWCzVjnXv3h2NRoPJZKqznoeHB0Cdz3rDDTegKAr33XcfR48e5aabbjprLH379qV9+/a8+uqrFBcX1zhf9XxWq7VGF0RwcDDh4eHVYvbw8GiyrgqtVlutJQTsY4PS0tKqHavr63auzy6Es0jLh2j1rrzySj7++GNuu+02Jk+ezB9//MGLL77I6tWrGThwIP/4xz+IjY0lLy+P7du3s2LFCvLy8hrs/sOHD+ef//wnL7zwAomJiYwdOxadTsfhw4dZvHgxb731Ftdcc805X69Xr17ccMMNvP/++xiNRgYPHszKlStrrMcBNPhzrlq1invvvZdrr72WTp06YbFY+N///udIdM4Us1ar5aWXXsJoNGIwGBg5ciTBwcGAvZVm/PjxLF68GF9f31qnDf+dRqNh3rx5TJgwgbi4OG699VYiIiJIS0tj9erVeHt78/PPP1NUVERkZCTXXHMNPXv2xNPTkxUrVrB161Zee+01x/X69u3L119/zZw5c+jfvz+enp5MmjSpXl+fc3X55Zfz9NNPc+uttzJ48GB2797NwoULadeuXbVy7du3x9fXlw8//BAvLy88PDwYOHAgMTEx5/TsQjiNE2faCNHkqqbabt26tca5V199VQXUyy+/XDWbzWpmZqY6c+ZMNSoqStXpdGpoaKg6atQo9aOPPnLUqZqCWtvU0ltuuUX18PCocfzJJ59Ua/uv99FHH6l9+/ZV3dzcVC8vL7V79+7qQw89pKanpzvKnMtUW1VV1bKyMnXWrFlqQECA6uHhoU6aNElNTU2tMaVTVdULes6/TzM9evSoetttt6nt27dXXV1dVX9/f/XSSy9VV6xYUa3e36faqqqqfvzxx2q7du1UrVZb67Tbb775RgXUO++886zPf7odO3aoV111lRoQEKAaDAY1OjpanTp1qrpy5UpVVVXVZDKpDz74oNqzZ0/Vy8tL9fDwUHv27Km+//771a5TXFysTps2TfX19XVM167ta6CqdX/vhw8frsbFxdU4Hh0drU6cONHxvry8XL3//vvVsLAw1c3NTR0yZIi6adOmWr//P/74oxobG6u6uLjUiONszy6Esyiq+re2PSGEuAj9+OOPTJkyhT///NMxxVcI0TxJ8iGEaBYuv/xy9u/fz5EjR5p08KcQouHJmA8hxEXtq6++YteuXfz666+89dZbkngI0QJIy4cQ4qKmKAqenp5cd911fPjhh7i4yN9MQjR38r9YCHFRk7+PhGh5ZJ0PIYQQQjQpST6EEEII0aQuum4Xm81Geno6Xl5eMrBMCCGEaCZUVaWoqIjw8HDHrsx1ueiSj/T09Fq3ORdCCCHExS81NfWsm0hedMmHl5cXYA/e29vbydEIIYQQ4lwUFhYSFRXl+D1+Jhdd8lHV1eLt7S3JhxBCCNHMnMuQCRlwKoQQQogmJcmHEEIIIZqUJB9CCCGEaFIX3ZgPIYQQLYeqqlgsFqxWq7NDEQ1Ap9Oh1Wov+DqSfAghhGgUFRUVnDx5ktLSUmeHIhqIoihERkbi6el5QdeR5EMIIUSDs9lsJCcno9VqCQ8PR6/Xy8KRzZyqqmRnZ3PixAk6dux4QS0gknwIIYRocBUVFdhsNqKionB3d3d2OKKBBAUFkZKSgtlsvqDkQwacCiGEaDRnW2ZbNC8N1Xol/yqEEEII0aQk+RBCCCFEk5LkQwghhKg0YsQIZs+e7ewwWjxJPoQQQohGsmvXLoYOHYqrqytRUVG8/PLLZ61z/PhxJk6ciLu7O8HBwTz44INYLJZqZdasWUOfPn0wGAx06NCB+fPnVzs/d+5cFEWp9urSpUtDPtoFkdkuQoiLQrHJwpebj9O7jS/92vo7OxwhLlhhYSFjx45l9OjRfPjhh+zevZvbbrsNX19f7rzzzlrrWK1WJk6cSGhoKBs3buTkyZNMnz4dnU7H888/D0BycjITJ07krrvuYuHChaxcuZI77riDsLAwxo0b57hWXFwcK1ascLx3cbl4fuVfPJEIIVq1OV8nsmxfJnoXDctmD6NtoIezQxINSFVVyszOWeXUTaet1ywNm83GQw89xLx589Dr9dx1113MnTu33vdduHAhFRUVfPrpp+j1euLi4khMTOT111+vM/lYtmwZ+/btY8WKFYSEhNCrVy+eeeYZHn74YebOnYter+fDDz8kJiaG1157DYCuXbuyfv163njjjWrJh4uLC6GhofWOuylI8iGEcLqk7GKW7csEoMJi4+ed6fxrVEcnRyUaUpnZSuwTS51y731Pj8Ndf+6/7hYsWMCcOXPYvHkzmzZtYsaMGQwZMoQxY8YwYcIE1q1bV2fd6Oho9u7dC8CmTZsYNmwYer3ecX7cuHG89NJL5Ofn4+fnV6P+pk2b6N69OyEhIdXq3H333ezdu5fevXuzadMmRo8eXa3euHHjaoxVOXz4MOHh4bi6uhIfH88LL7xAmzZtzvnr0Jgk+RBCON0fezKqvf/zcLYkH8JpevTowZNPPglAx44deffdd1m5ciVjxoxh3rx5lJWV1VlXp9M5Ps/IyCAmJqba+aqkIiMjo9bkIyMjo1ri8fc6ZypTWFhIWVkZbm5uDBw4kPnz59O5c2dOnjzJU089xdChQ9mzZw9eXl7n+qVoNJJ8CCGcbv3hHABmDG7L/I0p7EkrxGpT0WpkOe6Wwk2nZd/T485esJHuXR89evSo9j4sLIysrCwAIiIiGiyuxjRhwgTH5z169GDgwIFER0fzzTffcPvttzsxMjtJPoQQTlVhsbH9eD4ANwxowzfbUimtsHI0u5iOIc7/C000DEVR6tX14Uynt16APXabzQZQr26X0NBQMjMzq52vel/XWIzQ0FC2bNlyxjp1Xdfb2xs3N7dar+vr60unTp04cuRInbE3pebxL0EI0WLtTjNistjw99DTKcSTjsGe7DxhJCm7RJIPcdGpT7dLfHw8jz32GGaz2XF8+fLldO7cudYul6o6zz33HFlZWQQHBzvqeHt7Exsb6yjz22+/Vau3fPly4uPj64yruLiYpKQkbr755nN70EYm63wIIZxqS3IeAP3b+qEoCjGVs1ySc0qcGZYQtYqIiKBDhw51vqKjox1lp02bhl6v5/bbb2fv3r18/fXXvPXWW8yZM8dRZsmSJdXW3xg7diyxsbHcfPPN7Ny5k6VLl/L4448zc+ZMDAYDAHfddRdHjx7loYce4sCBA7z//vt88803/Pvf/3Zc54EHHmDt2rWkpKSwceNGrrzySrRaLTfccEMTfJXOTlo+hBBOtTWlKvmwr+0RE+gJwNHsYqfFJERD8PHxYdmyZcycOZO+ffsSGBjIE088UW2ardFo5ODBg473Wq2WX375hbvvvpv4+Hg8PDy45ZZbePrppx1lYmJi+PXXX/n3v//NW2+9RWRkJPPmzas2zfbEiRPccMMN5ObmEhQUxCWXXMJff/1FUFBQ0zz8WSiqqqrODuJ0hYWF+Pj4YDQa8fb2dnY4QohGZLWp9Hp6GUXlFn6+9xK6R/rw0850Zn25g37Rfnx792BnhyjOU3l5OcnJycTExODq6urscEQDOdP3tT6/v6XbRQjhNAcziigqt+Ch19I1zD6+o510uwjR4knyIYRwmi3JuQD0beuPi9b+46hqZdPckgqMZWanxSaEaDySfAghnGZrin2K7YC2p0b+expc8PewrwiZll/3rAIhRPMlyYcQwilsNpW/jtpbPgbEBFQ7F+Zj70s+aZTkQ4iWSJIPIYRT7Ek3kltSgafBhd5tfKudC/OxL5SUbix3QmRCiMYmyYcQwinWHswGYEiHAHTa6j+Kwn0rWz4KpOVDiJZIkg8hhFOsOGDfK2N4p+Aa56paPk5Ky4cQLZIkH0KIJncst4SdqQVoFBgTG1LjfFXLR7q0fAjRIknyIYRocj8mpgMwpEMgQV6GGufDfavGfEjyIURLJMmHEKJJqarKD4lpAFzRq/btyatmu2QYy7HZLqpFmEULN2LECGbPnu3sMFo8ST6EEE3qYGYRR7NL0LtoGBdXs8sFcLSGmK2qLDQmmrVdu3YxdOhQXF1diYqK4uWXXz5rnePHjzNx4kTc3d0JDg7mwQcfxGKxOM6fPHmSadOm0alTJzQaTbNMliT5EEI0qRX7MgEY2iEQL1ddrWUMLlp83e3nsopMTRabEA2psLCQsWPHEh0dTUJCAq+88gpz587lo48+qrOO1Wpl4sSJVFRUsHHjRhYsWMD8+fN54oknHGVMJhNBQUE8/vjj9OzZsykepcFJ8iGEaFLL99tnuYyuZaDp6YI87a0f2ZJ8tAyqChUlznnVc/9Um83GQw89hL+/P6GhocydO/e8HnnhwoVUVFTw6aefEhcXx/XXX8+sWbN4/fXX66yzbNky9u3bxxdffEGvXr2YMGECzzzzDO+99x4VFRUAtG3blrfeeovp06fj4+NzXrE5m0t9Cr/wwgt8//33HDhwADc3NwYPHsxLL71E586dHWXKy8u5//77+eqrrzCZTIwbN47333+fkJAz/6ARQlzcsotMbE3Jo39b/1oHiZ6LzMJydqYWADCqS80ptqcL9jZwOKuY7GKZbtsimEvh+XDn3PvRdNB7nHPxBQsWMGfOHDZv3symTZuYMWMGQ4YMYcyYMUyYMIF169bVWTc6Opq9e/cCsGnTJoYNG4Zer3ecHzduHC+99BL5+fn4+fnVqL9p0ya6d+9e7XfmuHHjuPvuu9m7dy+9e/c+5+e4mNUr+Vi7di0zZ86kf//+WCwWHn30UcaOHcu+ffvw8LB/Y//973/z66+/snjxYnx8fLj33nu56qqr2LBhQ6M8gBCi8eWVVDDx7XVkFZkI8jKw5J7BRPq5n0PFo7DxXTi2AcxlWF0iuVLTg2PhlxHsfeZt1qXlQzhLjx49ePLJJwHo2LEj7777LitXrmTMmDHMmzePsrK6Z2HpdKe6EjMyMoiJial2viqpyMjIqDX5yMjIqPHH+ul1Wop6JR9//PFHtffz588nODiYhIQEhg0bhtFo5JNPPmHRokWMHDkSgM8++4yuXbvy119/MWjQoIaLXAjRZP67Nskx9iK7yMQdC7bx3d2D8TCc4UdI4iL4eTZYTyUP4RzjDf0GTlZshKLF4FV3i2hV60pWoSQfLYLO3d4C4ax710OPHj2qvQ8LCyMry95dGBFR+wwtUT8XNObDaDQC4O/vD0BCQgJms5nRo0c7ynTp0oU2bdqwadOmWq9hMpkoLCys9hJCXDxOnxr71OQ4grwMHMgo4q4vEigxWWqvtP1/8MPd9sQjZjjc8BVlN/3KG9apFKpuhBXuhE/HQnFWnfetSj6yiyX5aBEUxd714YyXotQr1NNbL+yhK9hsNgAmTJiAp6dnna+4uDhHvdDQUDIzM6tdq+p9aGhorfc+nzrNUb1aPk5ns9mYPXs2Q4YMoVu3boC9SUiv1+Pr61utbEhISJ3NRS+88AJPPfXU+YYhhGhkhzKLySw04abTcv2AKLpF+HDTvM2sO5zDtHmb+frOQbjqtKcqHF0Lv8y2fx5/L4x5BjQa1u7J4C3zFLb4Xsoi1xdR8lPg65vg1t9Bo61x32Ave7eMdLuIi0l9ul3i4+N57LHHMJvNjuPLly+nc+fOtXa5VNV57rnnyMrKIjg42FHH29ub2NjYBnwS5zrvlo+ZM2eyZ88evvrqqwsK4JFHHsFoNDpeqampF3Q9IUTD2pxs3/a+X1s/DC5a+kb78eWdg/B117EztYC3Vh4+Vbg0D77/B9gs0O0aGPssaOw/Zlbst//11jWuF8pN34PBG1I3w6Z3a72vo+VDkg9xEYmIiKBDhw51vqKjox1lp02bhl6v5/bbb2fv3r18/fXXvPXWW8yZM8dRZsmSJXTp0sXxfuzYscTGxnLzzTezc+dOli5dyuOPP87MmTMxGE4N9E5MTCQxMZHi4mKys7NJTExk3759TfNFaADnlXzce++9/PLLL6xevZrIyEjH8dDQUCoqKigoKKhWPjMzs87mIoPBgLe3d7WXEOLisTfN3hXaK8rXcaxXlC8vXW3vF//4z6OkVe3B8tsDUJwJgZ3gincdzd1Wm8qqA1VTbIMhsAOMe95eZ9VzkJtU476OMR+SfIhmysfHh2XLlpGcnEzfvn25//77eeKJJ7jzzjsdZYxGIwcPHnS812q1/PLLL2i1WuLj47npppuYPn06Tz/9dLVr9+7dm969e5OQkMCiRYvo3bs3l112WZM924WqV7eLqqr861//YsmSJaxZs6bGKN6+ffui0+lYuXIlV199NQAHDx7k+PHjxMfHN1zUQogms++kPfmIC6/+h8G4uFDi2wWw6WguH/95lLlxmbDnO1C0cOWHoHNzlN1+PJ+8kgq8XV3o39Y+RozeN9nLH10Nq56Ba+dXu37VbBdjmRmTxYrBpWbXjBANbc2aNTWO/fDDD+d9vR49epxxau6MGTOYMWNGtWPR0dH89ttvZ7yuWs+1Sy429Wr5mDlzJl988QWLFi3Cy8uLjIwMMjIyHP1fPj4+3H777cyZM4fVq1eTkJDArbfeSnx8vMx0EaIZUlWVI1nFAHQK8apx/p5L2wPw3bYUrL8/Yj844E6I6FutXNWqpiO7BKPTVv7YURR7twwK7F0CGXuq1fF116HT2ltOcoorGuqRhBAXgXolHx988AFGo5ERI0YQFhbmeH399deOMm+88QaXX345V199NcOGDSM0NJTvv/++wQMXQjS+zEITZWYrWo1ClH/N6YqXdAikfZAHU6zL0OYcADd/GPFwtTI2m8ovu04CMCb2b92vod0g9gr75399UO2UoiiO1o+sQlloTIiWpF7Jh6qqtb5ObzJydXXlvffeIy8vj5KSEr7//vsWNT1IiNbkaI691SPKz+1Ui8VpFEXh5n7B3Odi/wNDvfQxcKs+in9LSh5pBWV4GVwY1bWWVU3jZ9o/7v6mxtRbGfchRMske7sIIeqUklMKQNvAupemnqqsJFAp5JgtmB3BV9Q4v3DzcQAu6x5WfUpulagBENEPrBWQsKDaqcDKlo9c6XYRokWR5EMIUafkypaPtgF1JB8WE+5b3wfgQ+skFm2tvoLl6oNZ/LzTfuzm+Oga1R36327/uPPLapuA+XvY98TIK5GWDyFaEkk+hBB1Sq5s+WgXVEfysfNLKEqnwj2U76zD+HlnOsZSMwCpeaXc/81OAKbHR9Mt4gy7b3adDDoPyEuC1C2Ow/6e9uQjt0RaPoRoSST5EELU6US+PfmobbApqgqb7K0eukv+RbtQf0wWG99uP8HO1AKufH8jeSUVxIV78+hlXc98I4MnxE62f75zkeNwgKPlQ5IPIVoSST6EEHXKqJxlEuZTyw60yWsh5yDoPVH6THd0qzzzyz6u/mAjOcUmuoR6Me+WfrWP9fi7ntfbP+77Eaz2PWP8PexjPiT5EKJlkeRDCFGrcrOVgsoulDBvt5oFtnxs/9jzBnD15uo+kUT62ctZbCojuwTz7d2DCfOppW5toi+xz5Qpy4fUv4BTLR8y4FSIlkWSDyFErTKM9lYPV50Gb7e/LYZccBwOVq7AOODOynJaFt0xiBmD2/LiVd2ZN70fnoZ6LKKsdYFOE+yfH/gVOH3AqSQfommMGDGC2bNnOzuMFk+SDyFEraq6XEK9XVH+viX51k9AtUG7ERDUyXG4TYA7cyfHcf2ANmg09dvGHIAuE+0fD/wCqlot+Wjuy0mL1mnXrl0MHToUV1dXoqKiePnll89a5/jx40ycOBF3d3eCg4N58MEHsVgsjvPff/89Y8aMISgoCG9vb+Lj41m6dGm1a3zwwQf06NHDsWdafHw8v//+e7UySUlJXHnllY7rTJ06lczMzIZ58LOQ5EMIUavMyuQjxPtv4z3MZbD9c/vnA+6kQbUfCS5u9paVzL0EVM52qbDaKDZZzlJZiItLYWEhY8eOJTo6moSEBF555RXmzp3LRx99VGcdq9XKxIkTqaioYOPGjSxYsID58+fzxBNPOMr8+eefjBkzht9++42EhAQuvfRSJk2axI4dOxxlIiMjefHFF0lISGDbtm2MHDmSK664gr179wJQUlLC2LFjURSFVatWsWHDBioqKpg0aRI2m63xviiV6rWxnBCi9ajqdqkx2HTPd1CWBz5toNP4hr2p3h3aX2rv0jnwC+4juuGq01ButpFXUoGXq65h7yeajKqqlFnKnHJvNxe3mq13Z2Cz2XjooYeYN28eer2eu+66i7lz59b7vgsXLqSiooJPP/0UvV5PXFwciYmJvP7669V2tj3dsmXL2LdvHytWrCAkJIRevXrxzDPP8PDDDzN37lz0ej1vvvlmtTrPP/88P/74Iz///DO9e/cGYNKkSdXKPPfcc3zwwQf89ddfxMXFsWHDBlJSUtixY4djN/kFCxbg5+fHqlWrGD16dL2ftz4k+RBC1OpkZfIRcnryoaqw+b/2z/vfDppG2Gm28wR78nF4OYz4DwEeBtIKysgtqSC6rsXOxEWvzFLGwEUDnXLvzdM2466rZbp4HRYsWMCcOXPYvHkzmzZtYsaMGQwZMoQxY8YwYcKEM+5SGx0d7Whd2LRpE8OGDUOv1zvOjxs3jpdeeon8/Hz8/Pxq1N+0aRPdu3cnJCSkWp27776bvXv3OpKL09lsNoqKivD39681JqvVyuLFiykpKXHsMG8ymVAUBYPB4Cjn6uqKRqNh/fr1knwIIZwj87QxHw4ntkHGLnBxhT7TG+fG7UfZP6Zvh9I8/D30pBWUkSczXkQT6dGjB08++SQAHTt25N1332XlypWMGTOGefPmOXZyr41Od6p1LiMjg5iYmGrnq5KKjIyMWpOPjIyMaonH3+vU5tVXX6W4uJipU6dWO757927i4+MpLy/H09OTJUuWEBsbC8CgQYPw8PDg4Ycf5vnnn0dVVf7zn/9gtVo5efJknc/XUCT5EELUKqO25CNhvv1j3FXgXvtfWRfMJwKCukL2fji6Bn+PSEBmvDR3bi5ubJ622Wn3ro8ePXpUex8WFkZWln3Tw4iIiAaLqyEsWrSIp556ih9//JHg4OobN3bu3JnExESMRiPffvstt9xyC2vXriU2NpagoCAWL17M3Xffzdtvv41Go+GGG26gT58+aDSNPxxUkg8hRK0y/97tUl4Ie+2719J3RuPevP1Ie/KRtJIAj38AssR6c6coSr26Ppzp9NYLsMdeNQizPt0uoaGhNWaPVL2va7f30NBQtmzZUu1YXXW++uor7rjjDhYvXlxrN4ler6dDhw4A9O3bl61bt/LWW2/x3//au07Hjh1LUlISOTk5uLi44OvrS2hoKO3atavz+RqKJB9CiBqsNtWxjb1jwOnuxWAuhcDO9p1oG1OHkfDXe3BkFf6d7gFkczlxcahPt0t8fDyPPfYYZrPZcXz58uV07ty51i6XqjrPPfccWVlZjpaM5cuX4+3t7egyAfjyyy+57bbb+Oqrr5g4ceI5xW6z2TCZav4/CgwMBGDVqlVkZWUxefLkc7rehZDkQwhRQ26xCYtNRaNAUOW29myv3O6+7y1Qj5kD5yV6iH1cSVE67TUn7DFJy4e4CNSn22XatGk89dRT3H777Tz88MPs2bOHt956izfeeMNRZsmSJTzyyCMcOHAAsLdGxMbGcvPNN/Pyyy+TkZHB448/zsyZMx2DQxctWsQtt9zCW2+9xcCBAx1jQdzc3PDxsW/g+MgjjzBhwgTatGlDUVERixYtYs2aNdXWA/nss8/o2rUrQUFBbNq0ifvuu49///vfdO7c+YK/Tmcj63wIIWqoGu8R6GnARauB9EQ4uRO0euhxfeMHoHOD6MEAdCnZBsiYD9H8+Pj4sGzZMpKTk+nbty/3338/TzzxRLVptkajkYMHDzrea7VafvnlF7RaLfHx8dx0001Mnz6dp59+2lHmo48+wmKxMHPmTMLCwhyv++67z1EmKyuL6dOn07lzZ0aNGsXWrVtZunQpY8aMcZQ5ePAgU6ZMoWvXrjz99NM89thjvPrqq438VbFT1Its2cDCwkJ8fHwwGo2OucdCiKa1bG8Gd/4vgR6RPvx07yXw24Ow5SPodjVc82nTBLH+TVjxJFnhoxhw9PZTsYhmoby8nOTkZGJiYnB1rWVjQtEsnen7Wp/f39LyIYSoodrqplYL7F1iP9HzhqYLou1QAPxztqLBJpvLCdGCSPIhhKjh5Omrmx5bDyXZ4OZv38ulqYT1BL0XLhWFdFWOSbeLEC2IJB9CiBoyTm/52P2t/WDsFaBtwuXNtS4QbV+NcZBmH2VmK2UV1qa7vxCi0UjyIYSooarbJcxTA/t/sh/sdnXTB1LZ9TJYa58JkCvTbYVoEST5EELUULWpXJeSBCg3gmeoY/ZJk2prH2A6QLMfDTbpehGihZDkQwhRQ2ahvYUhImu1/UDXyxtnE7mzCe0BBm+8KCVWSZG1PoRoIST5EEJUU1RupthkQcGG1/EV9oOdJzgnGK2Lo8UlXrNPNpcTooWQ5EMIUU3VeI+BhuNoijNB7+kYe3E+jCYjO7N3kpCZQIoxhXovLVTZ9TJIs1+6XYRoIWR5dSFENRlGe5fL5YZEqMC+yZuLoV7XsNqs/Jr8K98e+pYdWTuqnQt0C2Ry+8ncHHszgW6BZ79YZfLRX3OAbcWl9YpDCHFxkpYPIUQ1VdNsB6uJ9gOdxter/r7cfVz7y7U8tv4xR+IR4h5CW++2GLQGcspy+HTPp0xaMonFhxafvSUktAcmrSfeShluOXvq+zhC1MuIESOYPXu2s8No8ST5EEJUk1lYjjfFtDUfth9of+k51/3m4Dfc+NuNHM4/jJfei1m9Z7HimhWsuHYFP1/5Mxtv2Mibl75JXEAcxeZint70NI+sf4RyS3ndF9VoyQ7oB0BofsKFPJoQTW7Xrl0MHToUV1dXoqKiePnll89a5/jx40ycOBF3d3eCg4N58MEHsVgsjvNr1qxBUZQar6oN5poD6XYRQlRz0lhGvGY/GlQI7ATe4Weto6oq7ya+y0e7PgJgVJtRPBn/JH6u1bcN12v1jGozihGRI/hi/xe8kfAGvx79lcySTN4b9R7uOvdar18cOhCy1hBTsqPW80JcjAoLCxk7diyjR4/mww8/ZPfu3dx22234+vpW21zudFarlYkTJxIaGsrGjRs5efIk06dPR6fT8fzzz1cre/DgwWp7qAQHBzfq8zQkafkQQlSTYTQxRFPZvREz/JzqfLb3M0ficU+ve3hjxBs1Eo/TaTVabom7hY/GfISHzoNtmdu4d9W9lJprH9NhbTMEgK4Ve8Emq5w2R6qqYistdcqrvoOcbTYbDz30EP7+/oSGhjJ37tzzeuaFCxdSUVHBp59+SlxcHNdffz2zZs3i9ddfr7POsmXL2LdvH1988QW9evViwoQJPPPMM7z33ntUVFQfcB0cHExoaKjjpdE0n1/p0vIhhKgms7D8VPLR7uzJx09JP/FGwhsAPNDvAW6Ju+Wc7zUgbAD/HfNf/rn8n2zN2Mqj6x/l9RGvo1Gq/xB1jepFoeqGt1IKGbsgvPe5P5C4KKhlZRzs09cp9+68PQHFvfZWtdosWLCAOXPmsHnzZjZt2sSMGTMYMmQIY8aMYcKECaxbt67OutHR0ezduxeATZs2MWzYMPR6veP8uHHjeOmll8jPz8fPr2aCvmnTJrp3705ISEi1OnfffTd79+6ld+9T//Z79eqFyWSiW7duzJ07lyFDhpzzMzqbJB9CiGqsxjTaa06iKhqUtmfewn5/7n6e2vgUALfG3VqvxKNKz6CevD/qfe5Ydgcrj6/k3R3vMqvPrGplArzc2GLrwmjtDixH/8RFkg/RiHr06MGTTz4JQMeOHXn33XdZuXIlY8aMYd68eZSVldVZV6c7tf9RRkYGMTEx1c5XJRUZGRm1Jh8ZGRnVEo+/1wEICwvjww8/pF+/fphMJubNm8eIESPYvHkzffr0OY8nbnqSfAghHMxWG13LtoMOLCE90bnV3XVSVFHE/Wvvp8JWwfDI4czuO/u879snpA9PDX6KR9c/yse7P6ZbYDdGthnpOO/tqmOrGstodmBJWofLJfed972EcyhubnTe7pwBw4qbW73K9+jRo9r7sLAwsrKyAIiIiGiwuM5X586d6dy5s+P94MGDSUpK4o033uB///ufEyM7d5J8CCEcsopMDFDsm7i5tBt2xrIvbnmR1KJUwjzCeO6S52p0ldTXpPaTOJB3gM/3fc6TG5+kW2A3gt3tA+g0GoV9hh5gXYjuxF/2cR/OWO5dnDdFUerV9eFMp7degD12m80GUK9ul9DQUDIzM6udr3ofGhpaa/3Q0FC2bNlSrzoAAwYMYP369XWev9hI8iGEcMgwltNPcxAApXI7+9psTN/IT0k/oaDw8rCX8TH4NMj9Z/eZzdaMrezP28+j6x/lozEfOZKaHI/OFBrd8DYXQcZuCO/VIPcUoj7q0+0SHx/PY489htlsdhxfvnw5nTt3rrXLparOc889R1ZWlmP2yvLly/H29iY2NrbO+yYmJhIWFnY+j+QUzWdorBCi0eVln6S95qT9TdTAWsuUmkt5etPTAEzrOo1ewb0a7P46rY6Xhr2Em4sbm09uZsnhJY5zvp5ubLV1sb9JaT5/4YmWJSIigg4dOtT5io6OdpSdNm0aer2e22+/nb179/L111/z1ltvMWfOHEeZJUuW0KVLF8f7sWPHEhsby80338zOnTtZunQpjz/+ODNnzsRgsK80/Oabb/Ljjz9y5MgR9uzZw+zZs1m1ahUzZ85sui/EBZLkQwhxSupmADJ0bcDdv9Yin+75lLTiNMI8wpjVe1atZS5EjE8M9/a6F4DXEl4jpywHAH9PPZsl+RDNiI+PD8uWLSM5OZm+ffty//3388QTT1Rb48NoNHLw4EHHe61Wyy+//IJWqyU+Pp6bbrqJ6dOn8/TTTzvKVFRUcP/999O9e3eGDx/Ozp07WbFiBaNGjWrS57sQilrvXZ4aV2FhIT4+PhiNxmqLpwghGt+m/84k/uQXbA+cTJ97aw5cyyjJYNKSSZRby3l9xOuMiR7TKHFYbBZu/O1G9uXuY3zb8bwy/BWe+HEPiX+t4ifD/4HBBx5OlnEfF7Hy8nKSk5OJiYnB1dXV2eGIBnKm72t9fn9Ly4cQwiG4IBGAwqDap+u9s+Mdyq3l9Anuw+g2oxstDheNC3Pj56JRNPyR8gfbM7fj76Fnr9qWco07mIyQKfu8CNFcSfIhhLCzmIgqszf/WiIG1Dh9KP8QPyf9DNgXE1MUpVHD6RrQlas7Xg3Ay1tfxt/dBStaDrt2txeQrhchmi1JPoQQdumJ6DGTo3rjFda5xumPdn2EisqY6DF0D+reJCHd0+sePHQe7M3dy/GKDQAkauLsJyX5EKLZkuRDCAGAevwvABJsnQj1rb4oU1JBEstSlgHwzx7/bLKYAt0CuaP7HQAsO/kZKGY22braTx7bIPu8CNFMSfIhhADAnLIJgG22ToR4Vx9IVtXqMarNKDr712wVaUw3x95MmEcY+RXZ6Hw3s7ksCvSeUG6EzL1NGosQomFI8iGEAFVFc8K+quIhfSyuulOzSFKMKfyR8gfQtK0eVQxaA3f2sE9N1AeuIbe8HFubQZXBSdeLEM2RJB9CCMhNwqU8F5OqI8+7+iqKn+/7HJtqY3jkcLoGdHVKeFd0uIJwj3A0LsXo/DZTGibJhxDNmSQfQghItY/32Km2I9DXy3G4oLzAMcNlRtwMZ0QGgE6j45897a0u+oC1pPn1tJ84tgEq99wQQjQfknwIIeB41XiPzoSdNtj028PfUm4tp6t/V/qG9HVWdIB94zmNNRCNSzFfFx2qHPdRAFky7kOI5kaSDyEEHLcvq77N1omwysGmZpuZLw98CcBNsTc1+roeZ6PT6Ai2TgBgWdp3mKv2npGuF9GARowYwezZs50dRosnyYcQrV1JLuQeBuzTbKtaPlYcW0FWaRaBboGMbzvemRE6xBiGYrN4YTTn8HtA5fbiknyIi9iuXbsYOnQorq6uREVF8fLLL5+1zqxZs+jbty8Gg4FevXo1fpBOIMmHEK1d5WZyxzRRGPEkzMfe8lHV6jG181T0Wr3TwjtdoJcn5rwhAHxWehQbyLgPcdEqLCxk7NixREdHk5CQwCuvvMLcuXP56KOPzlr3tttu47rrrmuCKJ3DxdkBCCGcrHKw6VZrRwBCfVxJKkhiR9YOtIqWazpe48zoqgnw0FORPxCP4LUcKUljnZcvw4vyIWsfhHZzdnjiDFRVxVLhnCTRRa+pV7ehzWbjoYceYt68eej1eu666y7mzp1b7/suXLiQiooKPv30U/R6PXFxcSQmJvL6669X29n2795++20AsrOz2bVrV73v2xxI8iFEa1c53uMviz35CPNx5b2d3wMwLHIYQe5BTgvt7/w99GBzI1S5lBPqb3waGMzwogJ714skHxc1S4WNj+5b65R73/nWcHSGc98BecGCBcyZM4fNmzezadMmZsyYwZAhQxgzZgwTJkxg3bp1ddaNjo5m7177IOhNmzYxbNgw9PpTLYfjxo3jpZdeIj8/Hz8/v/N/qGZOkg8hWjOLCdJ3APbBpj5uOly0Nn5K+gmAazpdPK0eUJl8AO7lI9C5LWe7rZzdej3dU9bBoLucHJ1oKXr06MGTTz4JQMeOHXn33XdZuXIlY8aMYd68eZSVldVZV6fTOT7PyMggJiam2vmQkBDHOUk+hBCtU3oiWE1UGPxJKQ+li48rq1JXUWAqINg9mMHhg50dYTUBHgYAikrcmRA7gZ+SfmKRjxcvHNtoH/ehkWFsFysXvYY73xrutHvXR48ePaq9DwsLIysrC4CIiIgGi6s1k/+pQrRmleM9snx7AQphPq58d+g7AK7scCUumovr75NAL3vLR3aRiWldpwHwh4c7ORUFkL3fiZGJs1EUBZ1B65RXfaeJn956URW7rXJQ84QJE/D09KzzFRcX56gXGhpKZmZmtWtVvQ8NDT2fL2OLcXH9ZBFCNK3K8R5Jbt0B8PIqZPXJv1BQuLLjlc6MrFYhXvaZOPmlZjr4dKFnUE92Zu9ksZcnd6esh5C4s1xBiAtTn26X+Ph4HnvsMcxms+P48uXL6dy5c6vucgFp+RCi9VJVxzTb3ZouABg19pVOB4UNIsLz4mte9nXXodfaf2xlF5m4seuNAHzj5YU5+U9nhiZaiYiICDp06FDnKzo62lF22rRp6PV6br/9dvbu3cvXX3/NW2+9xZw5cxxllixZQpcuXard48iRIyQmJpKRkUFZWRmJiYkkJiZSUVHRZM/Z2KTlQ4jWKjcJSnPAxZVtFW2AAlJM9gW7JrWf5NzY6qAoCsHeBk7kl5FZaGJ09GiC9b5kUcCyrC1MlHEf4iLi4+PDsmXLmDlzJn379iUwMJAnnnii2jRbo9HIwYMHq9W74447WLv21Myg3r17A5CcnEzbtm2bJPbGJsmHEK1V5X4uhPfhRIEVjesJ8irScXNxY1SbUc6N7QyCvezJR3ZROTqNH1O73MC7uz5gkUFhYvYBCIk9+0WEqMOaNWtqHPvhhx/O+3o9evQ449TcGTNmMGPGjLPG0NLInwhCtFaVg01pM5AMYzk6H/uU2xFRI3DXuTsxsDMLqdx7JrPQBMDVXabiAuxyNXDowA/OC0wIcc4k+RCitaocbFoa2o9ikwkXb/tKipe3u9yZUZ3VqeSjHIBAt0Audbf3s3+XutxpcQkhzp0kH0K0RqdtJnfcvTtajyQ0LsX4GfyID493cnBnFuRlX+ujquUD4OpOVwPwszmbclORU+ISQpw7ST6EaI0qZ7kQ2JnkEj06b3uXy9i2Y9FpdGeo6HxVLR9ZReWOY/HdbiLcaqNIo7A88WNnhSaEOEeSfAjRGp023uNwdj4uXva9KC72LheAEG97y0fWaS0fGq2Oq9zsXS/fHv3JKXGJ2qmq6uwQRANqqO+nJB9CtEapW+0fowaRkL0JRVuBpzaInkE9nRvXOQiuXGgs87SWD4Apna5Bo6psr8jlqPGoM0ITp6laVKu0tNTJkYiGVLXWiFZ77hv11abeU23//PNPXnnlFRISEjh58iRLlixhypQpjvMzZsxgwYIF1eqMGzeOP/7444ICFUI0EKsFTibaP4/sT9KOl0ALPfyG1nsZameoavkoKDVTbrbiqrP/EAzpOoVhW55njbsb3+9ZwANDnnJmmK2eVqvF19fXsSeKu7t7s/j3Jepms9nIzs7G3d0dF5cLW6mj3rVLSkro2bMnt912G1dddVWtZcaPH89nn33meG8wGM4/QiFEw8reD+ZSMHhj8o3CqOwEYFSbMU4O7Nz4uOnQu2iosNjILjIR5V85Ldjdn2sM4awhnx+Tf+O+QY+j017c41dauqr9S6oSENH8aTQa2rRpc8GJZL2TjwkTJjBhwoQzljEYDK1+0xwhLlppCfaP4b1ZnboRNCZsZh9GxvRzblznSFEUQrwNpOaVkVFYfir5AIZ0mEzQ4U/Ippx1aesY2WakEyMViqIQFhZGcHAwZrPZ2eGIBqDX69E0wCrCjbLC6Zo1awgODsbPz4+RI0fy7LPPEhAQUGtZk8mEyXRq4FhhYWFjhCSEqFKVfET05YfDvwOgK+9JoKerE4Oqn0hfd1LzyjiRX0r/tv6O4y5xV3JZ4jss8PHml0PfSfJxkdBqtRc8RkC0LA0+4HT8+PF8/vnnrFy5kpdeeom1a9cyYcIErFZrreVfeOEFfHx8HK+oqKiGDkkIcbq07QCYw3uxLcu+l0s798HOjKjeIv3cADiR97fdRQPaM8nVviHemvQNGE3Gpg5NCHEOGjz5uP7665k8eTLdu3dnypQp/PLLL2zdurXOteofeeQRjEaj45WamtrQIQkhqpiKIWsfAH/pFEy2EmxmL/qF9HJuXPVU1dWSml9zJkXn2GvpbKrArFpZmrK0qUMTQpyDRp9q265dOwIDAzly5Eit5w0GA97e3tVeQohGcnInqDbwjmB51jYALEXdiA33dW5c9VTV8pH695YPgLgrmVRcAsBPh77jYEYRV7y3gU6P/U6/Z5dzxXsbWH84pynDFUL8TaMnHydOnCA3N5ewsLDGvpUQ4mwqx3tYw3uzOnUNUJV8NK+kv6rl40RBLWtI+LXlMp8uaFSVnXn7uGH+L+xMLaDCaiOnuIKdqQXM+GwLiakFTRu0EMKh3slHcXExiYmJJCYmApCcnExiYiLHjx+nuLiYBx98kL/++ouUlBRWrlzJFVdcQYcOHRg3blxDxy6EqK80e2vHroA2FJjyUa1ueCud6BDk6eTA6qeq5SO9oByL1VbjfFCPG4gvsy9CVuyymY7BnqyYM4xfZ13C6K4hWGwqL/1+oEljFkKcUu/kY9u2bfTu3ZvevXsDMGfOHHr37s0TTzyBVqtl165dTJ48mU6dOnH77bfTt29f1q1bJ2t9CHExqBxsulpjn2FmKe7M4HbBaDTNa/GnEC9X9FoNVpvKSWN5zQJxVzKpxH7c3WcLj1zWmQ7BXsSF+/DUFXG4aBQ2Hc1lp7R+COEU9Z5qO2LEiDOu7b50qQzwEuKiVJQJxlRAYW2hfUdbS1FXBg2ofRr8xUyjUWgb6M6hzGKSsourrfUBgLs/vXz742E7RIm+CG+fE4B97aEIXzcm9QxnyY40Fm0+Ts8o3yaPX4jWTvZ2EaK1SLe3ehwP7sTRwhRUVYOlpDMjOgU5ObDz076yqygpu6TW85v0oxhTYh8T8svRX6qdu76/fUr/r7tPUm6ufRkAIUTjkeRDiNbihH28x5oA++Bva2k7BsdE1Gw1aCZOJR/FtZ7/KLMTlxZZAFie8jtm66kVNvu39SfSz41ik4UV+zMbP1ghRDWSfAjRWlTOdFlqs4+FsBR15eHxXZwZ0QXpEFyZfGTVTD6Sc0o4lGsmq6wPgRYrhZZSNp3c5Div0ShM7GFPwpbvk+RDiKYmyYcQrYHNBunbMWoUdptOAjA6ekSzHu9wpm6XjUn2dTwOBU1gXGXXy29JP1crM7prCABrDmbXOmNGCNF4JPkQojXIOwrlRv5090ZVVFRTKHMvG+bsqC5IuyAPFAVyik3kFJuqnUtIyQfAs9NwxuMBwKrjqyiznFqUrE8bP/zcdRjLzGw7lt90gQshJPkQolWoXN/jJ3f7Jmy9AwcT4t18NpKrjYfBhZgAe2KxN736hpRVyUTfmEB6dr2GCLOFMlsFa0+sdZTRahQu7RIMwEoZ9yFEk5LkQ4jWIC0BM7Dd1T5Nflb8FKeG01C6RfgAsCft1AZyWUXlHM8rRVGgdxtflJ7XM77E3jXzx5GfqtWv6npZsT+riSIWQoAkH0K0ChXHtpLgaqBCa8PTxZe+oT2dHVKD6BZhXxb+9OSjqsulc4gX3q46CO7KhMqdbtelb6SooshRdmjHQHRaheSckjpnzQghGp4kH0K0dBYTmqw9rHO3L0k+Kno4GqVl/NevavnYdeJU8lHV5dKvrZ/jWKduN9CuwkyFamXl8ZWO416uOga1sy+yJl0vQjSdlvETSAhRp8O7NuGimvnTzT4+4pLIS5wcUcPpGemLTquQVlBGSo69a2XDEftMlwExp1ZuVbpfy4TKWS9/HPq+2jVGVo77+POQ7HQrRFOR5EOIFm7bhhVkaLWk6LVoFA3xYfHODqnBeBhc6Bttb+H483A2GcZyDmQUoSgwtEPgqYLeYUzw6wbAX9mJ5JXnOU4NqSy37VgeJousdipEU5DkQ4gWLDmnBH3WDja62We2dAvsho/Bx8lRNazhnewtF2sOZvPnoWzA3iLi56GvVi66503EmkxYUVmWfGoPqo7BngR66ik320g8XtBkcQvRmknyIUQL9vXWVHopSWyoTD6GhA9xckQNb1TXqm6TbN5eZd8w79LOwTULdrmcy8rsS6z/fuhbx2FFURzjPjYdzW3kaIUQIMmHEC2WxWpjWcIBojUn2eRmH2w6JKLlJR+dQrzo39YPi03lRH4ZOq3C1P6RNQu6ejMuzP782wsOkVGS4ThV1XVz+qwZIUTjcXF2AEKIxrHmYDbhpQfY46WnSKvBW+9Nt4BuTR6HJT+fiqQkTEePYs3Lx1ZchKqqaNzd0Xr7oI9ug75de3QR4SiKcl73eO7K7tzy6RbySip45opuhPm41VoutMc0+qzdzHZXV/5I/p0Z3W4FIC7c3hX198XKhBCNQ5IPIVqon3el01NJYmNlq8egsEFoNdpGv69qs1G6dRvFq1ZSsnEjpsNHzqmeS1AQ7oMG4XXpCDxHjkTjeu4rsHYK8WLDwyNRsa9cWqd2I7jsdzPbXV1ZevgHR/IRG25fL+SksZzcYhMBnoZzvrcQov4k+RCiBTJZrKzan8XrmiPMrxzvcUlE406xtRYUkLdoEcbvvsecllbtnC4iAn37drgEBaH19AKNBltZKdbcPCpSkjGlHMOSnU3hzz9T+PPPaDw88L5sAv4zZmBo3/6c7q85U9JRRe/OqNABPG/ax57Co6QVpxHhGYGnwYV2gR4czSlhb3ohwzoFnc+XQAhxjiT5EKIF2piUS5HJTHv3ZPYY7Lu/xoc3zhRba0EBOR9/TMGXX2Erta+lofH0xGvcWDyHDsN94ABc/PzOeA1beTlliYmUbNhI4a+/Yk5Pp2DxtxQs/hbPUaMInn0fho4dGyTewC5X0HfzDra6ubI8ZTkzus0AoGu4N0dzSth/UpIPIRqbJB9CtEBL92QQTi6H3EzYFC86+LQn1CO0Qe+hWizkf/U1Oe+8g9VoH6hp6NKFgNtuxWvs2Hp1m2hcXfEYNAiPQYMI+vdsyhISyF2wgOIVKyleuZLiNWvwnXotQbNmnTWROatO4xm3+mG2urmyNOlnR/LRPtC+CFtKbsmFXV8IcVYy20WIFujPQ9n00hw5NcW2gbtcTEePkjLtRjKffRar0YihY0ciP/yAmCXf4zN5cr0Sj79TNBrc+/cn6t13affrL3iNGQNWKwVffsXRSZMpWrHiwoL3CGCUbywaVWVPwSHSiu1dRG0rk4/kHEk+hGhsknwI0cJkGMtJN5bTS3PEsbjY4IjBDXJtVVXJ+/xzkq+8ivJdu9B4eRH65BPELPkerxEjznu2Sl0M7dsT+c7btPl8AfoO7bHm5HDi3n+R/vDDji6e8xHYZRL9yk0ALEtZBpxKPo7lnv91hRDnRpIPIVqYxFT7xmptPI+Q5eKCq+JC35C+F3xdW0kJaffNJvP5F1BNJjyGDKHdzz/hd8MNKC6N24PrMWAAMd9/T8Cdd4JGg/HHn0i57npMycnnd8HOlzG2cq+XZcm/AxATYE8+ThrLKauQZdaFaEySfAjRwuw4XoAeM+mGLAD6BfbAoL2wqaMVx4+Tcv0NFC1bBjodIf/3OFHzPkYX2rDjSM5Eo9cTPOffRC+YjzYoENPhw6RcO5XitWvrf7GA9oxyi7B3veTtJ604DV93Hd6u9iTqWJ50vQjRmCT5EKKF2XG8gG5KMhvd7HubXNJ2zAVdr2znTlKmXofp8GG0QYFEL1iA/403NngXy7ly79+fmO++w61vX2zFxaTeM5OCb789e8W/Cew0sVrXi6IoxFQNOpVxH0I0Kkk+hGhBzFYbu9IK6KE9wHZXe2vH4AtYUr14/QaO3Xob1oICXLt3J+bb73Dv07uhwj1vuuBgoud/hs+VV4LVysnH/4/s995DVdVzv0iXiae6XlL+ACDSzx2AtILyBo9ZCHGKJB9CtCAHThZRbrbh57Ufs6IQ4eJJW++253Wtwj/+IPXuu1FLS/EYMoTo+Z+hC6llwzYnUXQ6wp5/joC7/glAzjvvkv3aa+eegIT3YRSe9q6X3H2cKDpBuK99gO7JgrLGClsIgSQfQrQoO1LzAZV8j0wABgf1Oa/ukcLly0m7/wEwm/G+bAJRH7yPxsOjgaO9cIqiEDx7NiGPPw5A7rxPyH799XNLQDQaAjtNcHS9LD+2nHBf+1L06UZJPoRoTJJ8CNGC7DheQAQ5JBjsCceQDhPrfY2i1atJm3M/WK34XDGZ8FdeQdHrGzrUBuV/042E/F9lAvLxPLLfePPcKnaZyLjKrpelKUsdG9KlS7eLEI1Kkg8hWpAdx/OJNSSSotfhosKAyKH1ql/y11+kzbqvssXjMsKefx5F2/ib0TUE/xtvPNUC8tFH5M6ff/ZKMcMYZVbQqCp7c/eiM9inKadLt4sQjUqSDyFaiLySClJyS/H03ANAD70fXnqvc65ffugQJ+79F6rZjNeYMYS/9GKzSTyq+N90I0H3zwEg68WXMP7665kruBgIaDeS/pVdL/sL1wOQXWyiwmJr1FiFaM0k+RCihahaXMzoaV/f45Kgc5+VYs7MJPXOf2IrLsa9Xz/CX3sVRadrlDgbW8Add+B3880ApP/nEUr++uvMFTqfmvWy/uRK9C4aVBUyC6XrRYjGIsmHEC3EjuMFuFHCPlf76pyDO15xTvWsxcWk3vlPLBkZ6Nu1I/Ldd9Bc5GM8zkRRFEIe+Q9eE8aD2cyJ+2ZTcexY3RU6jmFUqcnR9RLiZ1/jQ7pehGg8knwI0ULsOF5Ae/ctlGg0+NtUurYZftY6qs1G+gMPYjp4EG1gIFEffYTW17fxg21kikZD+Isv4tazJzajkdR7ZmItLq69sLs/AVGDHLNe9D57AcgqMjVVuEK0OpJ8CNECWG0qiakFeFWO94jX+aPRnH28Rs5771O8Zg2KXk/UBx+gj4xo7FCbjMZgIOKdt3EJDqYiKYn0Bx9CtdUxjqPLRMZUdr2U6XcAknwI0Zgk+RCiBUjKLqbYZKHAMwOAIecw3qNo1Spy3nsPgNCnn8Kte7dGjdEZdMHBRL73LopeT/Hq1WS/807tBTtfxujSUhRVpUhNQnEpIKtIxnwI0Vgk+RCiBdhxPB9FW8gJg328R3yHyWcsbzp6lPQHHwLA76ab8J0ypbFDdBq37t0Je+ZpAHI/+JDidetrFvKLJtC/E30qu15cvPaQXSgtH0I0Fkk+hGgBdhwvINQzAYCuFWYCo+te38NWXk7afbOxlZTg1q8vIQ8/1FRhOo3PFVfge8P1AKQ//DDmzKyahTqMZkypvevFxXs32cWSfAjRWCT5EKIF2HG8AC+P3QAMcfEDl7pnq2S+9JJ9h9rAQCLfeKPZTqmtr5D//AdD165Y8/JIf+ABVIuleoGOYxldYp/h4uJ+jJNFmU6IUojWQZIPIZq5onIzh7KMGKv2cwnsWWfZwmXLKPjyKwDCX3wRl6CgJonxYqAxGIh4/TU07u6Ubt1KzvvvVy/QJp4QrRu9KrtestWtTohSiNZBkg8hmrldJ4wohpOUulhxt9no1W5creXM6emcfPz/AAi443Y8LxnSlGFeFAwxMYQ+9RQAOR98SGlCwqmTLnpoN8Ix66XCsFNWORWikUjyIUQzt/1YPh4e9im2A8vK0bUZXKOMarGQ9uBD2AoLce3enaBZs5o6zIuGz6TL8ZkyBVSV9P88grW45NTJDqMdyYfWPYWDOWnOCVKIFk6SDyGauR2pBXh42hfGGoI7eIXUKJM77xPKEhLQeHgQ8fprF/0utY0t5LFHcQkPw5yaStbLL5860XEMYVYr3ctNKIrK0uQVzgtSiBZMkg8hmrnd6ZmUumUDMDioV43z5QcPkV25nkfI/z2OPiqqKcO7KGm9vAh//gUACr75huK1a+0nfCIhONax18uGk6ucFaIQLZokH0I0Y9lFJvJt+1AVlWizmajoYdXOq2Yz6Y/8B8xmPEeOxOeKc9vvpTXwGDQQ/1umA5D++ONY8u0b89FhNKMrp9wmFe0ityzXWSEK0WJJ8iFEM7b/ZCEungcBGFxaDlEDqp3P+e9HmPbtR+vjQ9hTc1EUxRlhXrSC/v1v9O3bY83OIfOZZ+0HO44l0mKlk8mCio1VqdL6IURDk+RDiGZsX7oRV4/9AFxiViE4znGufN8+cj78EICQJ/6vVU2rPVcaV1fCX3wRNBoKf/uNojVroM0gTFoPJpTYN6JbnrLcuUEK0QJJ8iFEM5aQfgirvgidqtIvIA60LkBVd8ujYLHgNXYs3pdd5uRIL15u3bvhf8stAGQ89TTWsgoyAwc5xn1sydhCQXmBEyMUouWR5EOIZmxfgX0hrD7lJtyjBjmO586fj+ngQbS+voQ++YR0t5xF0L/uRRcZieXkSbLffJOiyEtpY7HQtkLBqlpZnbra2SEK0aJI8iFEM1VutpKv2tf3GFJWBlEDAahITSXnPfvqncEPP4xLQIDTYmwuNO7uhD41F4D8hQuxWSIBuLzYPgh16bGlzgpNiBZJkg8hmqm9J3PRuicBMKS0HCL7o6oqGXOfQi0vx33QIHymyOyWc+U5ZIhj8TH9R/PYb45ibOWsl83pmzGajM4NUIgWRJIPIZqpZUmbQGMmyGKho3c0uPtT+MuvlGzYgKLXEzb3Seluqafghx9C6++PLfkoRw+GEmO20F5xw6JaWJO6xtnhCdFiSPIhRDO1NXMTAEPKylGiBmEtKCDzBfvCWYF334W+bVsnRtc8ufj5EfzQgwBE7s/EXKJhbGEBAMuPyawXIRqKJB9CNFPHy3YA9uSDqP5kvvoq1rw89B3aE3D77U6OrvnyueIK3Pr2RW+1cGKHP+OM9kXGNqZvpKiiyMnRCdEySPIhRDN0sugkJk06GlUlvqyc0gIvjN9+B0DYU0+1+r1bLoSiKIT+3+PYFIXyE3pCUrW0c/HGbDNL14sQDUSSDyGaoV+O2Pci6WaqwEtxJfO9hQD4XHUV7n37OjO0FsG1Sxd29x8LQGaCN2NLKgDpehGioUjyIUQztOb4nwBcUlaGMbsd5fv2o/H0JHjOv50cWctx5LIbKDB4UFGkY9TqTAA2pG2gxFzi5MiEaP4k+RCimbHYLBwwJgAwpMBE9gb7MuBBs/6FS2CgM0NrUXyD/ZkXN8n+ZpcbPcq8qbBV8OeJP50bmBAtgCQfQjQzu3N2U6GW4mmFwARXrCUVGDp2wO+GG5wdWosS5GVgZVRfykJ9UK0abltRBkjXixANQZIPIZqZDWkbALj8eBnGI+4AhDz2OIpO58ywWpxATwMoCr8NngxA+J4y2qerrDuxjlJzqZOjE6J5k+RDiGbmz9T1oKqMX6OAquA1fhwegwY6O6wWJ8jLAMAvHgPw6Wg/9o9VGsotZaxLW+fEyIRo/iT5EKIZySvP40D+PobsU3HPcEHRKYQ89JCzw2qRAj3tyUduaQWB11+KorXRLtXMoAOqdL0IcYEk+RCiGdmUvgkXs42bV6sABI7viS483MlRtUwBnva1UsxWFVPPyQR0sc9yuWm1jU0paymzlDkzPCGaNUk+hGhGNqRt4LJtKv5FKi7uFvyn3+jskFosg4sWHzf7OJpM/34E9AAXNyvBRhjxVykb0zY6OUIhmi9JPoRoJmyqjV2H1nHlRhsAwT2K0LQb5OSoWraqcR9ZpaCJHU1Qj0IArtpgY+3un85YN7+kgqd+3svcn/aSXWRq9FiFaE4k+RCimTiYd5DRK3JxrwCDfwXu3QPBQ9b1aExBleM+sotN0PkyfNqWYQvW4F4BoYtWY7LWnlRYbSq3LdjKZxtSmL8xhbu+SMBmU5sydCEuapJ8CNFMbN/8E6MT7b/AQnsV4tKmv5MjavkCK1s+sotM0HEMilZHdLcsAIbtMLP5ryW11vvfphR2HC9wvE84ls+GpJxGj1eI5kKSDyGaCZ95P6BV4US0DvfgCpSIPs4OqcWr1vLh5gedxuEZXEFGZw+0KpS+93GNOjabyrz1yQA8fUUc0+OjAfgxMb3pAhfiIlfv5OPPP/9k0qRJhIeHoygKP/zwQ7XzqqryxBNPEBYWhpubG6NHj+bw4cMNFa8QrVL2nytpv68AiwaCetiXUydCNpBrbIFe9hkvOUX2jeXoNQ0Av55F2BSI3p5O4Y5t1ersSTdyIr8MT4ML1/aNYkxsCADrDmejqtL1IgScR/JRUlJCz549ee+992o9//LLL/P222/z4YcfsnnzZjw8PBg3bhzl5eUXHKwQrZFqtXLyhRcAWNXbjV4eudjQQFhPJ0fW8lVr+QDoMAbcA+ihyWBzT3tikvziM9WSio1JuQAMauePm15L/7b+uOo0ZBaaOJRZ3LQPIMRFqt7Jx4QJE3j22We58sora5xTVZU333yTxx9/nCuuuIIePXrw+eefk56eXqOFRAhxbow//IguOY1iV1jXKxIAk39n0Hs4ObKWr2q2S07VbBUXPfS8AQ1gHOJOhRb0Ow9Rsn69o84mR/IRAICrTsvAGPvn64/IuA8hoIHHfCQnJ5ORkcHo0aMdx3x8fBg4cCCbNm2qtY7JZKKwsLDaSwhhZzOZyH73HQCWDNbQ3qwAoGvTz5lhtRqBf2/5AOh/B6BwSUkSf/Szfz8yX30N1WbDYrWxLSUPgMHtT81E6hvtB8CeNGPTBC7ERa5Bk4+MjAwAQkJCqh0PCQlxnPu7F154AR8fH8crKiqqIUMSolnLX7gIy8kMcrxgRT89Y8vtMy1coiT5aArBlS0fucUmrFVTZf1joPME+pSbWDPYlRIDVBw8SOGvv3Isr5SSCituOi1dQr0c1+kW4Q1I8iFEFafPdnnkkUcwGo2OV2pqqrNDEuKiYC0qIve//wVg8VANAe7d6K/YZ1EQLjNdmoK/hx5FAZsK+aUVp04MvAstMNhSxA/x9h+j2W++RVKavdWjXZAHGo3iKN4t3AeApOxiSissTRa/EBerBk0+QkNDAcjMzKx2PDMz03Hu7wwGA97e3tVeQgjI/eQTrEYjOSFurO2uEFgWgbdShlnjCsFdnR1eq+Ci1eDvbh9YWm2V0phhENKNsUVGfu+nkO+lwZyWRvl33wLQIdiz2nWCvV0J9jJgU2H/yaImi1+Ii1WDJh8xMTGEhoaycuVKx7HCwkI2b95MfHx8Q95KiBbNkp1N3oLPAVhwiRmbRiE6x/4Xc2lAHGh1zgyvVaka95Fz+rgPRYGh99O33ISHRuWbS+yHI3/7GleLiQ5BnjWu0y3C3vqxN126XoSod/JRXFxMYmIiiYmJgH2QaWJiIsePH0dRFGbPns2zzz7LTz/9xO7du5k+fTrh4eFMmTKlgUMXouXK+eAD1LIyTF2i2dxRpa13DJ0KjwEy3qOpBZ2+yunpYq/AJbAzI0tKWdNdoSjYE7eSQiYfXU/74FqSj3B7q+7uE5J8CFHv5GPbtm307t2b3r17AzBnzhx69+7NE088AcBDDz3Ev/71L+6880769+9PcXExf/zxB66urg0buRAtVMWxY+R/sxiAtZPagKLQzW8QPZQkANxjBjgzvFbHMd22+G/Jh0YLwx5kbGkpVq3C15WtH9ccXkN7t5qLicU5Wj5kRp8Q9U4+RowYgaqqNV7z588HQFEUnn76aTIyMigvL2fFihV06tSpoeMWosXKfvsdsFjwGHoJ33seBCBEiaOrYm/5kGXVm1agZy1jPqp0u4r+ntH4WK0s71jGMR8/vMxleP+yuEbRziH22S9J2cWnZs4I0Uo5fbaLEOKU8n37KPz1VwCMt04mrzwPD50HvtkmDIqFUq03+MU4OcrWpc5uFwCNFt2YZxhZWoaqUfh+gH09j8L//Q9Lfn61olH+7uhdNJgsNk7klzZ63EJczCT5EOIikvX6GwB4T5zIn+7HARgcPhjdyV0A5Pl2tw92FE3m1IDTitoLdBzLGM92AOzumUpmSDS2khLyPvmkWjGtRqF95UDUw7LMumjlJPkQ4iJRsmWLfZluFxeC7pvFnyf+BGBoxFACCvYAoEqXS5M7Y8sHgKIwaPRLeFltFOtUDo+2T4PO+2IhluzsakU7Vg5EPZwlyYdo3ST5EOIioKoq2W+/DYDvtddgDHRjX+4+APqHDKa9+RAA3u0HOS3G1qrOAaen0UX2I7bCPqA0NXg9br16oZaXk/Pfj6qVO5V8yFofonWT5EOIi0Dppk2UbUtA0esJvOsu1qfZNyqLC4ijvMBCeyUdAO/2MtOlqVV1u+SVVmCx2uosl1IyHoANLsUEXGVPEgu+/hpzerqjTMcQe/JxRFo+RCsnyYcQTmZv9bBvHud7/XXoQkJYl7YOgGGRwyhI2opGUcnUBKN4Bjsz1FbJz12PVqOgqpBXUvu4j8JyM0dyeqKzashycSEp93+4DxiAajaT88GHjnIdgu0zXo5kFWOTGS+iFZPkQwgnK1m/nrLERBRXVwL/8Q/MVjMb0zcC9vEethPbAEhzlyXVnUGrUQjwsE+3zapj3EdSVjGoLmjKuwOw1JxN0PiOABR8/z0Vx+zTpKMD3NFpFUorrKQby5ogeiEuTpJ8COFEqqqS/ZZ9rIffDTfgEhTE9qztlJhL8Hf1Jy4wDs+cnQAY/Xs4M9RWLczXDYC0gtoThqTsEgBC9EMB+MPDA33653hcMhisVnLefx8AnVZDu0AZdCqEJB9COFHx6jWU79mD4uZGwB23AzhmuVwScQkaRUNosX3gqRre12lxtnaRfvbk40R+7clH1RiOXgED8TX4kOOiZYvVSNBQfwCMP/2M6cgRADpUjfuQ6baiFZPkQwgnUVWV7HfsYz38b7oRl4AA4FTyMSxyGBSm42/NwaoqeMZI8uEsp5KP2hcHS8q2JxKdQnwZ19Y+8PQXDw/c0r/Ca+RwUFWy33sPwLHpnMx4Ea2ZJB9COEnRihWY9u9H4+6O/223AZBamEpKYQouiguDwwdjOb4VgINqG9qEymBTZ4n0cwfqbvlIqmz5aB/kyeXtLgdghacHpRVFBA6xr3pa9PsflB886Nh0LjmnpLHDFuKiJcmHEE6g2mzkvPMuAH63TMfFz/4L6s80e6tH75DeeOm9KDm6CYBddCC4cr0J0fTO1O1SYbFxLM/eItIh2JOeQT2J9IykTIHV7m64nvgKr7GjAch+5x3aBXoAcDRbkg/ReknyIYQTFC1diunQITSengTMmOE4vjp1NQDDI4cDoFbOdDnhHodGI8uqO0ukb93dLsdyS7DaVDwNLoR4G1AUhcvb21s/fvELgrI8goYFgUZD8YqVhGelAJBbUoGx1NxkzyDExUSSDyGamGq1kv2uvf/ff8YMtD72lTELKwpJyEgA4NKoS8FqwTN3NwDGgJ7OCVYAEFHZ8lFUbqmRMBxxdLl4oFTuuzMxZiIAm3SQq9FgSF2M98TLACj+8H1CvO2tWEdzZNCpaJ0k+RCiiRX+9jsVSUlofHzwv2W64/j6E+uxqBba+7SnjXcbyNqHi7WcQtUNfWgXJ0Ys3PUujoThSHb1hKFqsGnVWA6Atj5t6R7YHSsqf/gFgfE4QRPiQKulZO2fXFKRAci4D9F6uTg7ACFaE9ViIady1kPArbei9fJynKvqchkRNcJ+IM3e5bLT1p62gZ40BzabSnmxGVOpGVOZhYpSC6YyC+ZyKzarDZsNVJtqX91TBa1Og4teg9bF/lGn1+LqqcPNS4+bpw4XvdbZj+TQKcSLzEIThzOL6Bvt5zh+5LTBpqeb2G4iu3N282tAKDfmZqI//i0+U67A+N33XLbtZ77rNl2SD9FqSfIhRBMy/vILFSkpaH198bvpJsdxs9Xs2M/l0jaX2g9WjvdIVDvQI8CjyWOtjaqqlBWZyc8oIT+jFGNWKcX5JvuroJzSgooGXTZcZ9Di5q3HO8AV70A3vANd8Q5wwzvIDf8wD3SGpktOOoV4se5wDof+tj5HVUtIh+Dqycf4tuN5Zesr7Dbnk6Iz0DZlHYHXzsH440+EH9lFXOhRjmaHNVn8QlxMJPkQoomoZjM579lXugy443a0nqcSiq2ZWyk2FxPgGkD3QPsS3erxv1CAHbYOTPJ3b/J4bVYb+RmlZB0rIvtYIdmpReRnlGIqtZy5ogIGNxf0lS+Dmwt6Vy0aFw2KoqDRKigaUFCwmG1YzVYsZhuWChtmk5Wy4grKi8zYbCpmkxVzdhmF2WVAfo37+AS6ERDpSUC4B4FRXoTEeOPh0zizgjqF1Fyfw2ZTScqyt178PfkIcAtgcPhg1qWt49foHsw8shV92q/4Xn01BV9/zfT9fzC/m4zlEa2TJB9CNBHjjz9iTk1FGxCA37Rp1c6tPn6qy0WjaKAoAyUvCZuqsJ0ujgGPjclssnIyqYC0QwWcPFJA9vEiLBW17OKqgHeAK36hHvgGu+MV4IqHrwFPP/vL3VuPRnthw8lUVaWizEJZkZnSQhOFOeUU5pRRmGv/WJBZSlmRGWN2GcbsMo7uyHbU9QpwJbSdD6HtfAhr70NgpCdKA8wU6hhi7yI7lHkq+UgrKKPMbEWnVWhTS4J4ebvLWZe2jl90Fu4BlF3fEHj7Wgq+/54euUfx2peIzTZUZjKJVkeSDyGagFpRQc77HwAQ8I870Lif+kWlqiprTqwBThvvccy+sdx+tQ3evoHoLvCXeW1sVhsZyYUc35NL2qECslIKa3SZ6Axagtp4ERTtRXC0F/5hnvgGuzX6WAxFUTC46zC46/ANcSe8Y80ypYUV5KYXk3uimNz0ErKPFZKbXkJRbjlFueUc3poJgKunjsgufkR19Seqqz9e/q7nFVOnEC80CmQWmjhpLCPMx409aUYAOod61fo9urTNpbi7uHOiPJcdgdH0yTmGLncTPlOnYly4kOv2/MZJ4x1E+DV9y5YQziTJhxBNoOD7JZjT09EGBeJ3/fXVzh3IO0BGSQauWlcGhQ2yHzy2AYAtti5EBzTcL6byEjPH9+aSsjuX43tza3ShePobiOjkR3hHX0Lb+eAb4n7R/lXu7q3H3dufqC7+jmMVZRYyUwrJOGok46iRk0eMlBebObItiyPbsgDwC3UnpmcgMb2CCIn2PudWEU+DC90ifNh1wsjmo3lM6R3Brsrko3uEb6113FzcGNt2LD8c+YEfwjvQJ+cYJMwn6J8Lyf7qG2LzjpG2bDUR1028sC+GEM2MJB9CNDJbRQU5H34IQOCd/0TjWv0v7zWpawCID4/H1aXyXGXLx1+2rhecfJSXmDmamM2RhCxOHMhHPa11w+DhQpvYAKK6+hPRyRfvwMbv3mlMejcXRwsHgNVqIzO5kNT9eZzYn0dmin3cSn7GcbYvPY6Hj56YnkG06xNERCe/syZag9oFsOuEkU1JuUzpHcHuE/bko0ekT511pnSYwg9HfmBp2Qn+o9XhfmILOjWHXX1H02/L72g++y/q1Msca4QI0RpI8iFEIyv4ZjGWjAxcQkLwnXptjfNVU2wvjaqc5VKaB1n2nWy32rpwt3/9Z7pUlFtI3pnDkW2ZHN+Xh816KuEIiPAgulsgbbsHENLO56Jt2WgIWq2G8A6+hHfwZeCkdphKzRzfl8fRxGyO7cmlxFjBnj/T2PNnGu4+ejr2D6HzgFACozxrTQbi2wXw0Z9HWX8kB6tNZdeJAgC6R9SdfPQJ7kMbrzYcLzrO8vYDuOLQBti+gMzLrqUsYSUeKYcpXr0ar5EjG+vLIMRFR5IPIRqRrbyc3P/+F4DAu/6JxlB9JkZ6cTr78/ajoNh3sQU4bt/P5bg2ijy8aXOOLR+qqpJxtJD9G9I5nJCFxWR1nAuI8KBDvxA69A3GN7j1ji8wuOvo2C+Ejv1CsJptnDiYz9HEbJJ2ZFFqrGDnilR2rkjFL8yDzgND6BIfVm32zMB2/ngaXEgrKOP15QcpLLfg5epC51CvOu+pKApXdLiCd3a8ww9ueq4A2PkVbYbfzk/tLuG6w6vIfvsdPEeMQNHIuo+idZDkQ4hGVPD111iys3EJD8P36qtrnF9xbAUAfUL6EOAWYD+YYl/vY7PVvqpp27Os8VFaWMHBvzLYvzGd/IxTe4/4BLvRsX8IHfuG4B9+cawTcjHR6jREdwsgulsAw67vxPG9uRzcnEnKrhzyT5bw1w9H2fJTMjG9AokbGkFkZz/c9S7cMCCKj9cl897qJABGdQk+64Dgye0n8+6Od9lWmESqXxRR+an0Kt/CMx2HMzllIxw4QNGy5XiPH9cUjy6E00nyIUQjsZWWkvPRxwAE3n03il5fo8zyY8sBGBM95tTBJHs3zOqKrgC1TuEEyDpWyM5VqRzZluXoVnHRa+jQJ5iul4QT1t5HxhGcI62LhpieQcT0DMJUZiFpexb7N6STcbSQpO3ZJG3PxifIjdih4dwVH8MfezNIzStDq1G4Y2i7s14/1COU+PB4NqZv5KfIrszMT6VN+h8U6W9gSfuhTDuwnOx338FrzGgU7cWzqqsQjUWSDyEaSf6XX2HNzUUXGYnvlCk1zmeWZJKYnQjA6Db2LdcpPAnZ+1FR2GiLI8TbgNtp01ptVhtJO7LZteoEGUeNjuPB0V7EXhJOx34h6N3kv/WFMLi5EDsknNgh4eScKGbvujQObs7AmF3Gpu+T2PpLMv/pE8LediqX9Aql2xnGe5zuivZX2JMPSw53A4bkFQQbruH7dsOYduIvKo4kUfjb7/hMurxxH1CIi4D8lBKiEViLS8idNw+AwHvuQdHpapRZeXwlAD2DehLiEWI/eHQNAPk+XSko92JI5aqZFeUW9q5LZ9eqVIrzTQBotAod+gbTY2QUIW29G/mJWqfASE+G39CZwVd14PC2THatPkHuiWJS/srEQ4GCbBvpqgthHc7eyjSyzUi8dF6kl+ewJSiGQdnJTPXazbumPhRMmorvl5+Q8957eE8Yj+IiP5pFyyb/woVoBPkLF2LNz0cfHY3P5Em1lllx3D7eo1qXy1F7l8sB934AdPRxZ+uvyexclYqpxL4mh5uXjrihEXQbHtFoS4mL6nQGLbFDwuk6OIy0g/kkrkjl2J5cknfmkLwzh5AYb/pd1pbobgF1JiGuLq5MiJnAN4e+4YfgCAZlJzOKrbxLH/YMGsfw37+jIiUF48+/4HvllKZ9QCGamAytFqKBWYuLyf30UwACZ95T61+xuWW5JGQmADCqzSj7QVV1tHxsMPViWJkLQWtz2fJzMqYSCz5Bblx6UxemPz+YgZPbSeLhBIqiENnFn8vv7ckNTw4k9pJwtC4aMpML+fW9XSx+YRtHE7OrraVyuikdpgCwovwkRYpCbNk2dFg4XGQj4I7bAch5/31Us7mpHkkIp5DkQ4gGlrdgATajEX27dnhPrH3lylWpq7CpNmIDYon0irQfzNpHWWEZ64v/gWF/WwaadGBRCYjwYOztcUx7ahCxl4TjopMBiRcD/zAPRzLYe0wbXAxaso8X8fuHu/n6ua0cSciqkYR0C+xGe5/2mGxm/vAPxWAtob/mAMk5JfhNm4Y2IABzaioFS5Y46amEaBqSfAjRgKxGI3nzFwAQdO/MOmcuVE2xrepyMZVZ2PztLv6X/SE7iy9Dqyqka21ccmtXrnt8AB37h7ToxcCaM3dvPYOv7sD05+LpOz4anauW3LRiln68h29e2Mrxvbmoqj0JURSFKzteCcD3fvZVWEdpdpCcU4LG3Z3AO/8BQM4HH2KrqHDOAwnRBCT5EKIB5c6fj62oCEPHjniNH19rGaPJyJaTWwC4NHQk25ce43+PbWTb7hDMqhu+vmUs9jDxa7CNHgNCZbpsM+HmqWfQlPZMf24w/Sa2ReeqJSe1mJ/f2cmPb+wgI9k+O2lS+0m4aFzYYy3igF7HKM12copNGMvM+F5/PS7BwVhOnqRg8WInP5EQjUeSDyEaiCU/n/wFnwMQ+K9761ytcnXqaqw2K8OKJ/Pnq+lsWpKEqdSCn0sq431fImBSMCk6Gx1Dal/iW1zcXD10DJzUjpufjafn6Cg0Lgpphwr47qUEfv9wN+QbHON8vvX2pq0mk/ZKur31w2Ag4K5/ApD734+wlZc781GEaDSSfAjRQPI+/RRbaSmG2K54jRlTZ7n1W3Zw9a4HiN09ilJjBV4Broy6tIDrA2bTvm0pe4rt02Y7VE6zFc2Tm6eeS67pyE1Px9NlcBiKAkcTs/nq6c0MODwZg9mdX728KFUUhml2kZxTDIDvNdfgEhaGJSuLgq+/dvJTCNE4JPkQogFYcnPJ+2IhAEH3/qvWFouCzFJ+fDeBmD+HE1gaiYubhiHXdODGpwbRhe/RKDbofBkHMgoB6BxS934hovnw8ndl1PSuXP9/A4npGYiqQvY2CzcmPkHbk5fwh7sXgzX7OJpdAoBGryfw7rsAyPnoY2ylpWe6vBDNkiQfQjSA3I/noZaV4dq9O56Xjqh2zlRqZv3iw3z59GZO7DFiw8qJ6N1Mf2YwvUa3QatWwJFV9sKdJ7D/ZBEAXcNk4bCWxD/cg8vu7sGUf/cmIMITvcWNS1KuJintRaLMLqRkFzrK+l55JbqoKKy5ueQvWuTEqIVoHJJ8CHGBzFlZ5H/5JQBBs061eqiqyoG/TrLwyb/YuTIVm1XFGJLG4p4v0XGSF26elXu9JK0Ecwl4R2D0iSWtoAyALqGSfLREEZ39mPpYf/pdG0GZSzHu5WGsLPgPUQkZFGTZWzkUnY7Ae+4BIHfeJ1iLS5wZshANTpIPIS5Q7kcfo5pMuPXqhccll9iPpRfzw+s7WDl/P2VFZvxC3bnkH9F81e4V8t0zGd/2tJkwe76zf4y7kv2Z9laPCF83fNxrLskuWgaNRmHgqM5kTF7PzrDVqFjwKPXgq6e3sPXXZKxmGz6TLkffti3WggLy//e5s0MWokFJ8iHEBTBnZDgGBQbdNwuzycrG747wzbNbST9cgIteQ/yV7bnu8QHsdN2AikqvoF6Ee4bbL1BRAgd/t3/e7SoOnLQ3vXcNk/EercGVcVewqe0P/Nz9BTz0h7FabGz5OZmvnt3CiSOFBN57LwC5n83HWlh4lqsJ0XxI8iHEBcj58ENUsxn3fv3IcG3Pl09tZsfy49hsKjE9A7nhyYH0GReN1kXD78n2JGNCzIRTFzj0B5hLwa8thPeR8R6tTHx4POFuQaR75uAV+SztJobj7q2nILOUn95MZHNmNGqn7tgKCx2L1wnREkjyIcR5qjiRRsF332PSe7M77g7++O8eivNNeAW4MvGeHlx2dw+8A9wASC1MZU/uHjSKhrFtx566yJ7v7R+7XQ2Kwv7KmS4y3qN10Cgaru58HQA/+bhi0B9k2lOD6D4iEhQ4vDWL9W3+SVrYJeQu+BxLfr6TIxaiYUjyIcR5yv7gA9ID+rIlfi7HjlnRaBT6jo/mhicH0rZHYLWyv6fYWz0Ghg4k0K3yXLkRDi+zf97taiosNg5k2Fs+YsMl+WgtpnS8Eo0K211dyU5bgcHNhWHXd+La//QjqI0XZovCwc43kNDhdo5/uNDZ4QrRICT5EOI85OxKYu2RcA50uRmzYiCojRfXPtqPQVPao9NX389FVVV+O/ob8Lcul/0/g7UCgrpAcCwHM4qosNjwcdPRNsC9KR9HOFGwezDdFfsYoATTtlPHo7255j/9uOTajrhoocCvE0uPdmD7Twew1bFrrhDNhSQfQtSDzaayc2Uq375/lHy/LmhUC/FXteeah/sSGFn7INF9uftIMiZh0BoYHT361ImdX9k/dr8GFIXEEwUA9IzylWXVW5lhwfbZTxsMhZSbyxzHNRqFnqOiuO7JQQSYT2DT6tn0WzpLXk0g76RMvxXNlyQfQpwjY3YpP7y2nfWLD2NFi2/BYa68MYg+Y6PRaOv+r/TDkR8AGNlmJF76ygQl/xikrAMU6HE9ADtTCwDoFenTiE8hLkaD46YRarZQqNWwbM//apz3DXZn4i1t6XxwIVpLORlHC/n6uS0k/JGC1WpzQsRCXBhJPoQ4C1VV2fNnGl89u5WTSUZcMNP50JcM9d9N6LDeZ6xbYa1wjPeY0n7KqRO7KvfsiBkKvlEAJByzDybs1ca3oR9BXOTahgQwuNA+OPnLQ7XvZus5dCgdQkoZuPUZQl3zsVlU/vrhKN+/nEB+hrSCiOZFkg8hzqA438Qv7+xk7aKDWExWQiN09P/raSJObiB49qyz1l97Yi1Gk5Fg92AGhg20H1RV2GlfEZWe0wA4kV9Kck4JWo1C/7b+jfU44iLlaXAhsCQWF1VlT3kGe3P31iijKApB983C1VRA7MqnuPTKcAzuLmQdK+Lr57aya3UqqowFEc2EJB9C1EJVVQ5tzeCrZzZzfF8eWp2GS67tSL+UBbiV5+E96XJcO3U663V+OvITAJPaTUKrqRyImroZ8o6CzgO6TgJgw5EcAHpF+eLlKiubtkYFnv0YU2JfXv2bg9/UWsZ90CDc+/cHsxm/DV9y/f8NJCrWH6vZxrqvD/PT24kU55c3ZdhCnBdJPoT4m/JiM8vm7WX5J/swlVoIjvZi6qP96eB5ktING8DFhaB//eus18kpy2Fd2joAJneYfOpEYuVGYbGTweAJwNpD2QAM6VB9iq5oPUqC+3B9YTEAvx39FaPJWKNMVesHQMF336ErSGfSv3oy7PpOuOg0nDiQz1fPbOHQlgxUVVpBxMVLkg8hTpO6P48vn9nMkYQsNBqF/pfHcNVDffELdSf79dcB8Jt6LfqoqLNe67ejv2FVrfQI7EE7n3b2g+Yy2PuD/fOeNwBQbLKw6kAWAKO7Bjf4M4nmITgoBPeyIDpWVFBuNfFz0s+1lnPv1w/P4cPBaiX7jTdRFIXuIyKZ+lh/gtt6Yyq1sPzTfSybt5fyYnMTP4UQ50aSDyEAq8XGxu+O8NNbiZQaK/ALdefqh/sy4PIYtFoNxavXULZzJ4qrKwF33XXW66mqyo9JPwIwuf1prR57vgeTEXzaQNuhACzbm0G52Ua7QA+6R8hMl9Yqyt+d7bZOjtaPrw9+XWfrRdD9c0CjoWjZMsoSEwHwC/Xg6gf7MGBSDIpG4UhCFl89u4W0g7Iqqrj4SPIhWr2CzFK+ezmBHcuPAxA3LIJrH+1PcLR9lVHVZiP7zTcB8L/5ZnTBZ2+d2JOzh0P5h9Br9IyPOW0H222f2D/2uxU0GlRVZf7GFACm9I6Q9T1asegAD3aoHZhYXIK7qpBSmMLmjM21lnXt1AmfKVMAyHz1VUeSotFq6D8xhmse7otviDslBSZ+eHMHf/2YJFNyxUVFkg/Raqmqyv6NJ/n6+a1kHy/C4OHChH92Z8S0ztVWKS389VdMhw6h8fIi4I7bz+naiyunS45rOw4fQ2VrRtp2SEsAjQ563wzAmkPZ7DphxE2n5caBbRr2AUWzEh3gznZbRzxUlUnF9oGnXx/4us7yQbP+hWIwULYtgeLVa6qdC472Zuqj/ek6JAxUSPj9GEte3U5hTlntFxOiiUnyIVolU6mZ5Z/sZdXn+7GYrER08uX6xwfQrndQtXKq2Uz22+8AEHD77Wh9zt4tUlhR6NjB9trO1546UdXqETcFPINQVZW3VhwG4Ob4aAI8DRf+YKLZCvI0kKaNoED14DpjAQCrU1eTUZJRa3ldaCj+0+1JbNbrr6FaLNXPG7SMvLkrY++IQ+/mQmZyIV8/u4XDWzMb9TmEOBeSfIhW52SSka+f3crhbVkoGoWBV7Rj8uzeePq51ihb8N13mFNT0QYGOn7Qn80vSb9Qbi2ng28HegX1sh8sy4fd39k/738HYG/1SEwtwOCi4Y6hMQ3xaKIZ02gUovw9SbR1oKPZTF+3cKyqle8Of1dnnYB//AOtjw8VR5Iw/vBDrWU69gvhusf6E9rOm4pyK8s+2cvKz/dTUW6ptbwQTUGSD9FqqDaV7UuPseS17RTlleMd6MpVD/ah34S2aDQ1x1rYysrIee99AALvuguN+9k3e1NV1dHlck2na06N4UiYD5YyCI6DqIGUVViZ+5N9IambB0UT7FUz8RGtTxt/e9cLwPU2+7+3bw99i9la+6wVrbc3AXfbB0Bnv/0OtrLau1W8A9248v4+9LusLYoCBzaeZPEL28g+XtQITyHE2UnyIVqF8mIzv76/i01LklBtKh37h3DdYwMIjam7GyXvf19gyc5GFxGB39Rr6yx3uoWJazlScATVpiNhTzvSCsqgohQ2vmsvMPhfWFWY/fUOjuWWEurtyn2jOzbEI4oWoE2AO9tV+7+HUZlJBLoFklOWwx8pf9RZx2/aNHQREViyssj7vOa+MFU0Wg0DJ7fjin/3xtPPQEFmKd++vI3da07ImiCiyUnyIVq8k0lGvn5uC8f25KJ10XDpTV0Yc1ssejeXOutY8vLI/egjAILum4Wi15/1Pql5pbyw4TN7/cIefJ+Qx8hX17DluzegNAd8o1G7Xc3TP+9l6d5M9FoNb9/QW1Y0FQ5t/N3ZaWuPDQVd/jGuj7GvgPu/ff+rM0HQ6PUEzb4PgNyPP8aSf+aptRGd/Lju8QHE9AzEZlH586tDLJu3l4oy6YYRTUeSD9FiqarKjuXH+eG17RTnm/AJduOa//Ql9pLws05pzfngQ2zFxbjGxuJ9+eXndL/Hf96I4rkTgEeH3sHAGH+wlNPmwMcAmOLv45nfj7Bg0zEAXr+uJwNiZB8XcUp0gDtFuHNcY1/EbqprJAatgf15+9metb3Oet4TJ2KI7YqtuJicDz44631cPXRMuKs7Q67pgKZyTZBvXthKzgnphhFNQ5IP0SKVl5j57YPdbPzuCDabSsd+wUx9tD+BkV5nrVtx7Bj5X9o3fgt+6EEUzdn/m+xJM7Ip+2cUxUoX3x5M73MJX905iM/jdhCq5HNS9afPT4F8uiEZgGeuiOPyHuEX9pCixWnjbx/nsc3aAQC/jL1c3s6e/P5vX91dKopGQ8gDDwCQv+hLTEeTz3ovRVHoNboNVz7QB08/A8asMr59KYF969OlG0Y0Okk+RIuTmVzIN89tJWVXDhoXheE3dGLM7XHoXevuZjld1utvgMWCx7CheAwadE51PtlwCJ2ffUGof/S8BQClJIeBqZ8C8IFmGiVWF7xdXXjnht7cHN+2/g8mWrxIP3cUBbZY2tsPnNjGzbH2WVarjq8itSi1zroegwfjOWIEWCxkvfzyOd8ztJ0P1z02gOhuAVjNNlZ/cYCV8/djNlkv5FGEOCNJPkSLoaoqO1em8v2rCY7ZLNc81I9uwyPPeeXQssREipYuBY2G4Mq/JM+mxGThj5Tf0biUEGAIYWSbkfYTq5+DiiII68Uj/5nLz/dewqZHRjGpp7R4iNq56rSEers6ZryQvp32XtEMDh+Misqi/YvOWD/44YdAp6N4zRqK12849/t66ph4Tw/ir2yPolE4uDmDxS9sJTe9+EIeR4g6SfIhWoSKMgtLP9rD+sWHsVlV2vcOYupjAwhqc/ZuliqqqpL58isA+Fw5BddOnc6p3or9mSg+9t1rb+l2Iy4aF8jcC9sX2AuMfwE3g47ukT54GM6t9UW0XlH+7iSp4ZhdvMBcCpl7HK0fS44sobii7oTAEBOD/403ApD54gs1Fh47E0Wj0GdcNFP+3Qt3Hz35GaV8++I2Dvx18sIeSIhaSPIhmr289BIWv7iNpB3ZaLQKQ6/ryLg7u2E4w2yW2hSvXEnZ9u0orq4EzZp1zvW+3LUKrWsGLoqBqzpeBaoKSx8F1QaxV0D04Po+kmjFov3dUdGQ7hlnP3BiK0PCh9DOpx0l5hKWHFlyxvqB99yN1s+PiiNJ5H9V9/LsdQnv6Md1jw0gsosflgobK+fvZ/UXB7CYpRtGNBxJPkSzdiQhi8UvbaMgsxRPPwNXPtCHHpdG1XuDNtVsJuvV1wDwn3ELupCQc6pns6nsLbFvfX5p+ET7Pi6Hl8PRNaDVw+in6hWHEFWDTg/qutgPnNiKoijcFHsTAAv3L8Riq7tFQ+vtTdB99uQ55513sBYU1DsGd289k2b1YsCkGFBg3/p0+94wubI3jGgYknyIZslmtbHh28Ms/XiPfW+Wzn5MfbT/GRcNO5OCb7+lIiUFrb8/AXfccc711iTvAvf9qKrCzL63gNUMyx6znxx0N/jLsumiftoE2JOPrRb7jBdStwBwebvL8TP4kVacxrKUZWe8hu8112Do2BGr0Uh25Sq99aXRKPSfGMOke3ti8HAh61gR3zy/leP7cs/rekKcrsGTj7lz56IoSrVXly5dGvo2ohUrLazgp7cSSVxhH/nfe2wbJs/qiZvX2RcCq421sNCxeVzgzHvQenqec91PdtsXFfO29aa9Xzv7Muo5h8A9AIbef17xiNatquVjdbF9rQ/yk6EkBzcXN6Z1nQbAJ3s+OeN0WMXFhZBHH7FXX7QIU1LS+ccTF8DUR/oT1MYLU4mFn9/ZybbfUlBtMh1XnL9GafmIi4vj5MmTjtf69esb4zaiFco4auSb57eSdqgAnUHL+Du7MfiqDmi05/9POee997Hm56Nv3x6/qVPPuV56cTq7jWsAGBxwLZQbYc0L9pMjHgHX82uFEa1bdIAHAEeKdNgCKwc9n9gKwA1dbsDdxZ1D+YdYl7bujNfxiI/Hc+RIsFrJfP6FC1q7wzvQjase7EPsJeGgwuafjvLbh7sxlda+54wQZ9MoyYeLiwuhoaGOV2BgYGPcRrQiqqqyZ+0Jlry2nZICE36h7lzzn3607xN8Qdc1HU0mb+FCAEL+8x8U3bkvdb5g7wJUrFiKOzC6fV9Y/yaU5kJgZ+h76wXFJVovP3cdnpWzoooCe9sPVna9+Bh8mNrZniB/svuTs14r5OGHUPR6SjZsoGjpmbtqzsZFp+XSm7pw6c1d0LpoSNmVwzcvbCPnhEzHFfXXKMnH4cOHCQ8Pp127dtx4440cP368zrImk4nCwsJqLyFOZ6mwsmrBftZ+ecgxjfaa//TDP8zjgq+d+dKLYLHgOWIEnkMvOed6eeV5fHf4ewAqckfQJ9ACm/9rPzl6LmhlSq04P4qi0D7I/m/7uPupGS9Vbo69GZ1Gx/as7WzPrHvJdQB9dLRjDFPmCy9gKym54Phih4Rz1YN98PJ3pTC7jO9e2sZBmY4r6qnBk4+BAwcyf/58/vjjDz744AOSk5MZOnQoRUW17xnwwgsv4OPj43hFRUU1dEiiGSvMKeO7VxI48FcGigLxV7Zn3J3dznm10jMp/vNPStb+CTqdfXGmeli4fyEmaznWskiCdXEE7/ovmEsgvDd0nnDBsYnWrWOIfX2aHWpVt8s2sJgACHYPZnL7yQDM2z3vrNcKuPMf6CIjsWRmkv3++Q0+/bvgaG+mPtqfNnH+WMw2Vszfz9ovD2K12Brk+qLla/DkY8KECVx77bX06NGDcePG8dtvv1FQUMA333xTa/lHHnkEo9HoeKWm1r18sGhdju/N5Zvnt5KTWoyrp47J9/Wiz7joek+jrY1qNpP54ksA+N90E4aYc5+VUlRRxJcH7Hu/VOSOYHg4sKXyl8CIR6EB4hOtW6cQ+6Dnv4qCwD0QLGX2BKTSbd1uQ6NoWJe2joN5B894LY2rKyGP22dg5S34HNPhww0So6unjokze9J/YlsA9qxNY8lr2ynOL2+Q64uWrdGn2vr6+tKpUyeOHDlS63mDwYC3t3e1l2jdVJvKtt+S+fndnZhKLQS3tf+VFdml4XaAzf/ySyqOHkXr70/gPXfXq+4X+7+gqKIIN8KxFMVyo+V7+y+HiH78f3v3Hd5U9QZw/JukSdp07z0oo2zKKoIgeygg4kARFMXxw4UDBVzgBkRFxYGKiIrKcqDiANl779ICBVqge7dpm3XP749ApVKgxaYDz+d58qS987y9TfLmnkXT/jVWRum/q9nZOx9HMo3Q6Dr7whPry9dHeEQwIHIAAJ8d+Oyyx3Pv1Qu3fn3BaiX9lVdrbOI4tVpF3NBoBj/SFr3ByT6v0hs7OJ2QWyPHl65eDk8+iouLSUpKIjg42NGnkq4CphILv805wLafT4CAlj1CuHmCvX65plhzc8n64EMA/J94HI171YdgLzQX8vUh++yilux++FNAq9Sl9pW9n5V3PaQacS75OJltxBp5LvlYV2Gb+9vY23KsOLmCo3mXv5sR9OyzqJydKdmxg8JffqnR8ka18eO2ZzvjG+ZGaZGFn9/by+4VyXJ2XOmiajz5ePrpp1m3bh0nT55k8+bNDB8+HI1Gw8iRI2v6VNJVJudMMUum7eTk/mw0Tmp639Wc3qOao9HW7L9p5ttvoxQWom/RAq9bbqnWvgviF1BkKSLSPZq8rOY8rP0Vtc0EYXHQuG+NllP67wr2dMZd74RVEaR4drIvPL0DzH83GI3xiaF/ZH8Ego/3fXzZY2pDQ/F7yH6XL+PNmdhquHG/p78Lt07sSPNrghACtvyQxJ+fHsRcVvX5ZaT/jhpPPk6fPs3IkSOJiYlhxIgR+Pr6snXrVvz9/Wv6VNJV5MiOdJbO2ElBViluPnr7mALX1vzsryW7d1Pwvb2XStCUF1FpNFXet8BUwNfx9rseXX1H4k8hd2pW2Vf2mizvekg1RqVS0eRsu49DpT7gGQGKFU5WHDPpoXYPoULFyuSVl237AeB77z3ooqOxZWeTOfOtGi+3k05DnzEt6DmyGWqNiqQ9WSydvpPctH/fy0a6utR48rFw4UJSU1MxmUycPn2ahQsX0rhx45o+jXSVsNkUNiw6wsrP47GaFcJb2IdJD4is+bY/wmol/eVXAPC89RYM7dtXa/+v4r+i2FJMU++mmPNb8YDTcvSYIawzNO5T4+WV/tuaBdirXo5mFv/dlijx9wrbNPVuysCogQBVuvuh0ukIfsU+31D+kiUYt22vwRKfPYdKReueYQyf0AFXL3357LhJuzNr/FxSwyXndpHqjLHAxLJZe9i/5jQAHQdFMuSxWFzcrmyY9MvJ++YbTImJaDw9CZhQvaHP88vy+eawfTCyh9o9xJHjJxmt+cu+suckeddDqnHNguzJR2JG0d/dt4/8aZ81+Tzn7n6sSlnF4ZzDlz2uoVMnvO64HYC0KS+ilDmmd0pQtCcjnutMaDMvLCYbf3x6kM3fH0Oxye64kkw+pDqSdiyfxW/sIO1YATpnDdePa8M1NzVGrXbMh7glI4Os994HwP/pCTh5e1dr/08PfIrRYiTGO4ZOftfRM2cRBpUJS1AsNOnngBJL/3XnutsezSiGqB6gdYWiVEjbV2G7aK9orm9kT04+2le1cTwCJkzAKSAAS3IK2R9+WLMFP4/BQ8eNj8cS2z8CgD0rU/j5/b2UFJoddk6pYZDJh1SrhBDsX3OKn97ZQ0mBGe9gV257tjPRsY5tE5Q5YwZKSQku7dpVu5Hp6aLT5eN6PNnxSQ4cPc5dGvtQ1do+clwPyTFizvV4yTFSKrTQuLd9xT+qXgDGtRuHWqVm7am1HMw+eNlja9zdCXppKgA5876gLD6+xsr9T2qNmmtvacLAB1qj1Ws4k5jPkmk7SD9R4LBzSvWfTD6kWmMx2/jri3g2LDqKogiadAzg1kkd8Qo0OPS8xRs3Ufjb76BWE/TSVFTq6v3bz94zG6ti5Zrga+gW0g3Nto9xVZk47RIDTQc4qNTSf52/ux4/Nx2KgIT0Qmg+2L7i0A8XVL008mzEkOghALyz650qdXF179MH90GDwGYj7YUXEVbH9kpp0jGAWyd3wivQQHGeiR/f3s3B9Wdkd9z/KJl8SLUiP7OE72fs4sj2DFRqFdfe2oQB97eqkWHSL0UxGkmfMgUA79GjcG7Rolr7x+fE89uJ3wD7XQ9VaR4d0u2j9Z5q86i86yE5jEqlolWIfWbkg6mF0HwIaPSQfQTSD1yw/aOxj6JT69iRvuOyM96eE/T8c6g9PSmLjyfniy9qtPyV8Ql25bbJnYhu749iFaz7NpHVXx3GarY5/NxS/SKTD8nhTu7PZsm0neScKcbFXcuwJ2KJ7RdRI8OkX07mu+9hSU1FGxJCwOOPV2tfIQTv7HwHgBsa3UBL35aUbpiNiyjlkBJJVLdbHVFkSSrXOtTe6+vQmQJw9oBm9p4tHFhywbbBbsHc2eJOAGbtmoVNufwHupO/P4ET7fMaZb8/u8aGXr8UnYsTgx5sTdfhjVGpIGFLOt/P3EVhdqnDzy3VHzL5kBxGKIJtvxxn+Uf7MZdaCWzkwYjn4ghtVr3GnleqZM8e8hYsACDolVdQu1ZvFtyVySvZlr4NnVrHY+0fg6IMtNvt3Rm/dxtFsJdjq4skqXX5nY+z7SPa3GZ/PrAUKkku7m9zP+46d47lH+OX41UbxdTz5uG49rwOYbGQOmkywmKpkbJfikqlosPASG58PBZnNy3Zp4pZPG0HKYdyHH5uqX6QyYfkEGVGC79+uJ+dy08C0LpnKMMndMDNW18r51fMZtJeeBGEwHP4cNy6X1ut/UssJczcOROAsW3GEuYeBmvfwMlWyh6lCc5tbnREsSWpgtah9uTjSHoxZqtib2Pk4m3v9XJ05QXbe+o9eaDNAwB8sOcDyqyX70arUqkIfvVVNGerX7I/nlOzQVxCWHMf+7g+UR6YjFZ++WAfO387gVBkO5CrnUw+pBqXdaqIJWe/xWi0avre04KeI2PQONXev1vOnDmYk5LQ+PkROGlitfefe2Au6cZ0QlxDGNt6LGQmIHZ/BcBrllH0aRFY00WWpAuEebvg4eyE2aZwNLMItM4QO8q+cue8Sve5s8WdBLkGkVGSUT42zeVoAwIImmpvG5X9ySeUHriwTYmjuPs4c/OEDrTqEQICtv18gt/mHMBU4vg7MFLdkcmHVKMStqbx/Zu7KMwuw8PPmVue6Ujza2p3UsGyhASyP7XP9Bn0wgtovLyqtX9yYTLzD80HYGLcRFycXGDlFFRC4XdbZ5JcWtM+onaqjqT/NpVKVX7349CZs3OxdBprfz66AvJTLthHr9HzaOyjgD2Jzi7NrtK5PG64AY8brgebjdRJkx02+FhlNFo1vUY1p8/dzdE4qSu0E5OuTjL5kGqEzaqw7rtEVs0/jM2iENHKh9ue7Yx/RNVnjK0JitlM6sRJYLXi3r8f7gOr1xVWCMGrW1/Foli4NvRa+oT3gePr4Oif2FQaZljvoGczfzQOGgxNkv7pXPJR3u7DtzFE9wIE7Jpf6T5DoofQ0rclxZZi3t/9fpXPFfjii2j8/TAfP07WrFn/ruBXoEW3EG5+xj6LdUFWKUun7+TI9vRaL4fkeDL5kP614jwTP72zm4PrzgDQ6YYoBj/SDmdXba2XJeu99zAdOYLGx4egl16qdo+aH47+wLa0bThrnHku7jlUQsCKFwD4VXs9J0UwfZoHOKLoklSpViH2Hi/7T583KFen++zPO+eBqeiCfTRqDc/GPQvAj8d+5EBW1apRnLy9CXntNQByv/yK4vXr/0XJr0xApAcjnutMeEsfrBaFlfPi2bDoCDY5LPtVRSYf0r+SejSPxdN2kH68EJ2LE4MfbkuXG6MdNkz6pZTs2EHuPPtYBcGvvoKTr2+19k83pvPWTvtMn4+2f5QIjwg4sBjS96Po3HmpcAhqFfRsJmdolmpPh7NVfIdSCyiznO3h0nww+DaB0ryLtv2IDYhlaPRQAKZvn44iqvbh7dazJ9533QVA6uRnsWTW/oRwzm5ahjzajk43RAGwf81plr2zB2OBqdbLIjmGTD6kKyKEYM+KFH6atZfSQjO+oa7c9mwnotr61Ul5bMXFpE6abO/dcustuPftW639hRC8tvU1ii3FtPVry+gWo8FSCqteBWBnxFjy8KBTlA9eBsdMfCdJlQnzdiHQQ4/FJth7Kt++UK2B7k/Zf948G0yVt414suOTGJwM7M/ezy9JVet6CxDw9AT0zZtjy80lddIkhFL7dx3UahVdbozmhofaoHPWkJZUwOLXd5B6LL/WyyLVPJl8SNVmKrHw+5wDbP7hGEIRNO0cyC0TO+EVUHfjXmS8Mc0+mFhYGIGTn632/suSlrHu9Dqc1E683O1lNGoNbP4ACk+DZzjvFvUB4PrWQTVddEm6JJVKRacoHwB2nsz9e0XbEeAVCcYs2PRupfv6G/wZ124cYB94rNhctQacar2e0HfeQeXiQsmWreR8NvdfxfBvNGrnz23PdsYnxJWSQjPL3tnDvtWn5LDsDZxMPqRqyUwuZPEbOzixLxu1k4qed8bQf2xLtHpNnZWpcMUKCn74AVQqQmZMR+NWvcHEUgpTeGPbGwA8EvsITbybQMFp2PA2AAXXPs+WFCMAg2TyIdWBzpH2qpcdJ/P+XqjRwgB7+ww2vQ95Jyvdd3SL0UR5RJFTlsOHe6s+g60+uhFBL9jbO2W9/z4le/ZcUdlrgleggVsndaJp50AURbBx8VFWzovHYpLDsjdUMvmQqkQIwcH1Z84Og/x3N9rW14XWyjDpF2M+fZq05+1vkL7334+hY8dq7W9RLEzeMJlSaymdAjtxb6t77StWvAjWUoi8lp+t1yAEtI/wItjTpaZDkKTLOnfnY3dyHrbzB+BqMRQaXQc2E/w8HiqpHtFqtOWNT785/E2VG5+CffRTj8GDwWbjzFMTsOblXX4nB9HqNfQf25LuI5qiVqs4uiODpTN2kpdurLMySVdOJh/SZZnLrKycF8+6bxNRrIKotn7c9mxnAiI96rRcwmzmzFMTUIqKcImNxX/8Y9U+xsd7P+ZA9gHcde5M6zHNXt1ycqN95lCVGgZN58e9qQAMblO745VI0jktgj1w1ztRZLKy99R5CYBKBUPeBScXOLEOdlRePdIttBtDoocgELy05SUsStUG8FKpVAS9/BK6yEisaWmkTpiAsNXd3QaVSkW7PuEMe6o9Bg8dualGFk/bSeI22R23oZHJh3RJualGlk7fydEd9tlou97cmBsealMn3Wj/KfOdWZTt34/a05PQt99Cpa1emXak72DuAfub9ZSuUwhyDQKbFX6fZN+g4z2c1DZmd0o+ahXc2C6kpkOQpCrRqFX0a2kfVffns8lwOd/G0P8V+88rp0Dm4UqP8UznZ/DSe3Ek7whfHvqy6ud2cyN09vuoXFwwbt5C1vuzryiGmhTSxIsRz3cmNMYLq8nGX1/Es/rrw1jk7LgNhkw+pItK3JbOkuk7yEsvwdVTx01PtqfDgMg6rWY5p2j1GnLnzwcg5I3X0YaGVmv/zJJMnln3DALBsMbDGBQ1yL5i93zIOAjOXtD7BX7cYx+75NomfgR4ONdcAJJUTeeS3+UH0jBZbeSXmFmbmMkfB9NJixkF0b3tVYULR0Fp/gX7+zj78EznZwCYs28OKYUXjo56Mc7NmhH8qr3nV84nn1C0atW/D+hfcvXUc+Pj7ek8OApUcHhTGkun7yQ3TVbDNAQy+ZAuYDHbWLMggb++iMdqVghr7s2I5+MIaepV10UDwJKaSuqz9jpsnzF3V7tbrUWx8PS6p8kpy6GJVxOe6/KcfUVRRnnXWno/jzD48NNee/Jxc4fqJTeSVNO6N/UjyMOZ7GIznV/7i/avruSeL3YwbsEuus1Yy3TDBIRHGOQmwff3g+3CqpWW7r0J1rXBZDPxyIrnUKrRhdZzyGC87z47/sekyZhOnKix2K6UWq0ibmg0wx6PxeVsNcySaTtI2JpW10WTLkMmH1IFOWeKWTJtJ/EbU0EFnQZHMXR8LAaP+jG2hVJWxunxj6MUFODcujUBEyZU+xjv7HyHPZl7cNO68W7vdzFoDSAELH8KyvIhqC10Gsuu5DySc0ow6DQMbCV7uUh1S6tRM3VoS1QqKCyzIgRE+7vSMtgDIWDOzkJm+kxFODnDsZXw4//s1YhAdrGJR77dTf9Z6zkaPwihaDlp3Me931evCiXwmWdw6dgRpbiY0489hq3owtFV60JYcx9uf74zYc29sZoVVs0/zKqvDsveMPWYStSzztKFhYV4enpSUFCAh0fdNmj8LxFCcGhDKhuXHMVmUTB46Og3tiXhzX3qumjlhBCkPfscBT/9hMbLi6ilS9GFVe+OxO8nfmfievsst+/2fpe+EWfvmhz6EZbcA2oneHAtBLXhiYV7+GlvKrd2DOOt29rVbDCSdIUS04tIzjESG+FFgLu9KnBlfAYPLdiFVRHMbJvGbccmg2JBRHRlfevXmLAin+xi++ig/VoEYHFdx27jlwhFx+S2nzG6Y4cqn9+SmcnJW2/DmpmJ67XXEv7JHFROTg6JtboURbDr95Ps+PUEQoB3sCsDH2iFb4hbXRftP6E6n98y+ZAoM1pYuyCBpD1ZAES08qXvmBb15m7HOblfLyDj9ddBrSbi87m4du1arf33Zu7lvj/vw6yYGdt6LE92fPLsgU/Apz2hrAB6ToLez5FVZKLb9FVYbIJfHu1OmzBPB0QkSTVn8Y5TTPx+PwDPRCZxf9Y09EoJZULLXNsN/OVzJ6/f0ZVWIZ4oQqH/dyPJtMTjZG7M5jFLcNFVvcF26cFDJN91F6K0FO877yRoyouOCuuKnEnMY8XnhygpNOOkVdPjjma06BZcL9qrXc2q8/ktq13+49KO5bPo9e0k7clCrVFx7a1NGPJI23qXeBi3bydj+nQAAp55ptqJx6miUzy+5nHMiple4b0Y3368fYWlFBbfZU88wuKgx9MAfLstBYtN0D7CSyYeUoMwonM4z9/QAoCZyY3pX/o6W5UWOKssPOq0jB9tj9Iq9QewWVGr1MwZ+CYoeqy6JCavqvrgYwAurVsR8uYMUKnI+/Zbcr9e4IiQrlhojDe3vxBHeAtvrBaFNV8n8OdnBykzVq2LseR48s7Hf5SiCHb/cZLtv55EKAJPfxcG3N+qzsfuqIz59GlOjrgdW24uHkOHEvLmjGp9gykwFXDX73dxouAELXxaMH/QfHs7D8UGS8dC/E9g8IP/rQfPUApKLPR4czWFZVbeuyOWYbGysanUcBw8U8Cv+9Mos9hoFuDGMJe9uK572d4QFcC/hX1k1Kb9ePK3OfyV9SEIJ5YOXUyMb9NqnStn7lwy33ob1GrCP/4It549HRDRlROKYPeKZLb/fAJFEbh56+l3T0tCY7zrumhXJVntIl1SUW4Zq+bHc+ZIPgDNugTSc2QMOuf6UW97PltBASdH3on5+HH0LVsQ9c03qF2qPsroqfx8Hlr5MMklBzCofYjVTKFFQDiju0Tgt+FF2P4pqLVw14/QqAcAb/6RwEdrk4gJdOe3x3ugqYMZeiWpRlnN9tlv1023z4QL0LgPpb1epsvyVxEuCQQ7N2H5bYvRqqte/SKEIO2FFyj4/gdUBgORX87HpU0bBwVx5TKTC1nx+SEKMktBBR0GRBA3NBqNk7z5X5Nk8iFVSgjBke0ZrP8uEXOZDSe9hp4jm9H8mvo5cqcwm0l54EFKtm3DKSiIqEWL0AYGVHn/X/efYvLGCahcDyNsekqS/4diCkGFwqvO3zGa5fYNb/kc2twKwNGMIobM3ojJqvDZ3Z3of3ZgJ0m6KpTmwfq3YNsnoFhApWZl5CieVLah0pRxT8uxTOj8ZLUOKcxmTo0bh3HzFjReXkR++w366GgHBXDlzGVWNi05Svwmezdc/wh3BtzXCq/AupsQ82ojkw/pAmXFFtZ+m0jS7kwAAht50O/elnU6E+2lnN+zRe3qSuS33+AcE1Pl/Xen5DD658fRuO8DoaWd09NEubUh1N2JJlufo79lNQAZ3V8jsJ99WPYyi42bPtxEQnoRPZr68dXYONlATbo65R6HlVPh8M8AfOwaw0cBpQB80u8TuoV2q9bhbMVGUu69l7IDB3AKDibqu2/RBtXP7ulJuzNZsyABU4kVJ52aHiOa0eJa2Ri1JsjkQ6ogJT6HVV8epqTAjFqtovOQRnQYGIFaU39vOWZ98CHZH3wAGg3hc+bg1qN7lfcts1i5bt4jlDpvRoWG2X3ep2f4dVCSC0vvheNrsaHmGfODrDf057sHuhDuY+CpxXv57UA6fm46fnu8R3k3Rkm6au1fDL88ARYjT/uE8aenGm+9Dz8M+x4/F79qHcqam0vyqNGYT5xA17gxkQu+xsm7fratKM4r46/5hzmTaK+CatTOj16jmte7hvYNjUw+JMA+UumWH5I4sPY0AN5BBvrd27JeNio9X+5XX5Pxhn2K+6CXXsL7jturvK8QgpHfP8ch468gVLzSdRrDYwbDiQ2w7GHITwGtgeKhn3Lrak8S0ovQadQY9BrySyxoNSq+HBtHt8bVe+OVpAYr8zBiwc2YitK4LSSEkzoNXYO7Mqf/HNSq6n1BsaSmcnLknVgzMnBu25aIeZ+jcaufY2wIRbDnrxS2LTuOYhM4u2npOTKGJh2rXrUrVSS72kqkHy9g8es7yhOPNr3DuO25up+J9nLyf/qpPPHwe+zRaicer26ebk88gFsin2R4aBz8OA6+HGJPPLwbwf1/4dZ2KAvu70L3Jn6YbQr5JRYCPfTMGd1RJh7Sf0tAC1T3rcTi1ph3MzPQK4ItaVuYd3BetQ+lDQkh4vO5aDw9Kdu/n1P33Y+tuNgBhf73VGoVHQZEctuznfANdaOs2MKfnx1kxdyDlBXLLrmOJu98XGUsZhvbfj7OvlWnQIDBU0ffu1sQ0cq3rot2WUV//cXpx58Amw2fMXcTMHlylethhRBM2z6N7xK+A8C39A7WdApB9ddL9iHTUUGnsdB3Crh4Vdhvz6l8FEXQOtQTZ62mpsOSpAZBlOSS8M4NHHZOYaq/LypUfNj3Q3qE9aj2scri40m+dyxKQQEusbGEz52Lxs3VAaWuGTarws7fTrLrj2SEInDx0NF7VAyN2vnXddEaFFnt8h+VejSP1V8lUJBlbzjW/Jogrr2tKc6u1Ztqvi4Ub9jA6YcfQVgseA4fTvDrr6FSV+3GnCIU3tj2BosSF4FQ4ZHem19d9+Gds8e+QVAbGPIuhHVyXACSdBX4a98JWHoPGwNS+N7DDTe1nm9uXEy0Z/V7r5QeOkTK2PvsCUj79oR/9lm9TkAAMk4Wsmp+PHnpJQDEXBNEjxFN0Rvq/3tofSCTj/8Yc5mVrT8dL69icfXS02tUDFFtGkb1QdHatZx5bDzCYsG9f39CZ71T5bkibIqNV7e+yvdHvwdUXJveiNmlG9FiBZ0b9H4e4h4ETf0bw0SS6htFEQx4exUPFL7Lb2GJ7HZ2JlLryTe3LMdTX/2RfksPHSLl3rEohYW4dOhA+Cdz0Li7O6DkNcdqsbH9lxPsXZmCEODqqeO6kTFEx8q7IJcjk4//kFMJuaxdkEBhdhkALa8NptutTdG7NIwP26JVqzj9xJNgseA+YAChb7+FSlu1bxkmm4lJ6yexKmUVatQ8lCUYV5xsX9l0IAyZBZ5ydFJJqo6vtyYz5af9PO/xLYsC95Lu5EQHvT9zbl2Oi1PVB/g7p/TgIVLG2hMQfYsWRHz2KU5+9f+LUfrxAv6aH28fmAyIjvWnx+3NcPPW13HJ6i+ZfPwHlBaZ2bT0GInb0gFw89HTZ3QLwlvWn1loL6dwxQrOPDUBrFbcrx9E6JtvVjnxKDQXMn71eHZl7EKr0vB6Vj7XF+VRiBtuN72Nut3tIPvtS1K1lZptdJu+irwSM590WM3LxX9SpFHT0zWCWTf/VK0RUM8pO3yYlPsfwJaTgy4ykoh5n6MNrf9fDKxmGzt/O8meFSkoikDnrKHr8Ma06hGKSo58fAGZfFzFhCI4vCWNzd8fw1RiBRW0uS6Ua4Y3rpfDo19M3uLFpL/0MigKHoMHEzJjepWrWrJKshj31ziO5B3BTaXl/dTTdC4zsVtpQkrfj7mpZ5yDSy9JV7e3VyQye/Uxov1cmdb2dx458ysmtZoh3q15feg31e6CC2A+eZKUsfdhSU3FKTCQiM/nom/SxAGlr3nZp4tZsyCBzJOFAARFe9JrdAy+IfWzG3FdkcnHVSo31cjabxNIO1YAgG+YG71HNSewUcP5OwkhyP7oI7JnfwCA5y03E/zKK6g0VetlciTvCI+teoxUYyp+aJhz5jQxZgtfWAcy1+VeVk8agN5J9liRpH+joMTCgHfXkVFo4vpWgQxzns1k4y5sKhUBxiZ0DJjC8ze0xrOaDTEtGRmk3Hcf5mNJaDw9CfvwAwydGkZDcEURHFx3mq0/HcdisqHWqGg/IIKO10eh1cn3HJDJx1XHYrKx6/eT7FmZgmITOOnUdLkxmra9w+r1KKX/JGw20l99lfyFiwDwHfc//B9/vMrdadedWsfE9RMpsZYQYRPMSU0lTKXndfU45hZ04vkbWvDAdfVvTglJaojWH8nivi93YLEJVCjc5z2TpYG52FQqtIVN8TA9wtJx3fFzq14bCGteHqf+N46y/ftRabUEv/E6nkOHOiiKmleUW8b6hUc4uT8bADdvPdfe2pTGHfz/80O0y+TjKiGE4NjOTDb/cIziPBMAUW39uO6OZrj7NKyhv21FRZx5+mmM69aDSkXgC8/jM2pUlfYVQvBV/Fe8vfNtBIK4MjPvZGTi6RXFhg6zuOtXI54uWjZN7oObvuFUPUlSfbc5KZsXfzpIUpaRGF8nHnJ5kZfdSrCqVIiiZkSKR1j84HW4VvN1p5SWkjpxEkUrVwLg98gj+D36SIP58BZCcHxvFpuWHKMo197YP6y5Nz1GNMMnpH53J3YkmXxcBbJSitiw+Eh5FYu7jzPdb2tKo1i/BvMCPcd04gSnH3kU8/HjqPR6QmbMwGPQwCrtezg9l6mbXuVw8V8A3FJYzPM5uWibDcI6bA6DPz1AYkYR4/s04akBVZ94TpKkqlMUgVqtgsI0Ns7vwxPuakxqNbaSCLq5Pc1no3tV+31JKApZ77xDztzPAfAYOpTgV19B7dxwvlhZzDZ2/5nMnj9TsFkV1GoVbXqH0XlIowbT47AmyeSjASstMrP15+PEb0wFAU46NR0HRRLbLwKnBlivWLxxE2eeegqlsBCnwEDCPvgAlzatq7Tve2u382nCy6hdUlAJeCY3j9GFRai6jYd+L/HR+hO8+UciXgYtayb0wttVTgolSQ53Zhc7vh3G476eFGnUKGYfbg9/iSmDel/R4fIWLyb95VfAZkPfsgVh772HLjy8hgvtWIXZpWxccpQT++xVMS7uWuKGRtPy2uAGVTX+b8nkowGyWmwcWHOGXX+ctPdiAZp2DqTr8MYNrooF7O07sj+eQ/ZHH4Gi4NKuHaGz30cbULVJm95YvYxvTkxD7WTEVYE3MzPpWmJhqrifdkMfJb/EwvQ/EhACZtzShts7Rzg4IkmSyu1fwolfxvFwoD+ntVqEzZn/tZjCY12vrO2GcetWzjz5FLa8PNSenoS+OQO3nj1ruNCOl3Iohw2Lj5KfYR8h1TvIQNfhjYlq2/DuWF8JmXw0IIoiOLItnW0/Hy9v1+EX7kaPEc0IaepVt4W7QpbMTFKfmUjJtm2AvUdL0NSpqHWXvzOhCIXJq2bx2+kvUakEMWaFdzPSCHXyYLrH83ySElJh+zFdI3npxlb/iRe2JNUrq14hb9MsxgcGsNdZhxAqbm96L891HY9GXf27tJa0NE4//gRl+/eDSoXfQw/h98jDVe4JV1/YrAqHNpxhx68nKTPaJ6gLaepFt5ubNKieiVdCJh8NgBCC5IM5bP0piZwzRsDeajpuaDQx1wTZ61cboOKNm0idOBFbbi4qg4HgqVPwHDasSvumFafx1JrJHMzdDcDNRSU8l5ON3qcp3LkIi1cjPl1/nG+2JuOkUTOuZ2NGxoXLxEOS6oKiwKJRmBN/4yW/EH5xt7dx6Bx4DW/1moGPc/UHPFTMZjKmTSP/u4UAuHTqSMj0GejC6v+AZP9kKrWy+89k9q06hc2iANC4QwBxQxvhE3x1NkqVyUc9l3o0n20/Hyf1aD4AeoMTHQZF0rZXWINs1wGglJSQ+fY75H3zDQD6mBhCZ81CH93osvsKIVh+Yjmvb32dYksxOgWm5OQwrNgIzQbB8Dng4u3oECRJqq6yAvi0F+Qe5z335nzmY0GlthBgCODVbq/SLbTbFR22YNky0l9+BaWkBLWbG0FTXsRj6NAG+UWjKLeM7T8fJ2FbOghABc06B9J5cCO8Ag11XbwaJZOPeir1aD7bfz3BmcQ8ADROatr2DqPDoMgGMfPsxZTs2EHqc89jOXUKAK+RdxA4aVKVWq1nl2Yzbds0ViSvAKBtmYlpWTmEq3SoBr4BHe+Rw6RLUn2WcQjm9gNLCa+oB7AoKAe1PguAEc1GMKHTBAza6n/ImlNSSJ04idK9ewHwuOEGAl98ASfvhvlFJPt0Edt/OVHeKFWlgpguQXQaHIWn/9WRhMjko575Z9Kh1qho0S2YjtdHNcjGpOcoRiOZ771H3tcLQAicgoMJfu1V3K699vL7CoUfj/7A2ztmUmQtQSME4/ILuD+/kJLoG/AY8gb4XP6uiSRJ9cCBpfD9fQCMtYwnpVkxGdi7x4e6hfLiNS9ybejl3xf+SVitZH/yCdkffQw2GxpvbwKfexaPIUMa5F0QsA+jsP3XE+WDlKnUKmK6BNJ+QGSDr46RyUc9IITgVHwuu1ckcyYxHzibdFwbQsdBkQ066RBCUPTHH2RMn4E1IwMAr9tuJWDixCpNl308L4mX1zzB7qKTALQwmZmanUNeWXO8Bz1Hq243OLL4kiQ5wu+TYNscioQLt5tf5PbbYvjm+JukGdMAGBQ1iImdJ+JvqP7U9KX79pH2wguYjh4DwLV7d4JemoouLKxGQ6hNGScL2fHrCZIP5pQvi2rrR4eBkQQ39qzDkl05mXzUIZtN4diODPasPEXOmWLg6kk6AEzHT5Dx2qsYN28BQBsWRtDUKbj16HHZfY0WI59te5Mvk37ACrgoCg/lF6PLiWWBuIGnR99E3xaBDo5AkiSHsFngq2GQvIlc4cbLvjN4fdxNfLj3Q75N+BZFKLhp3Xi8w+Pc1uy2aveIEWYzOZ9/TvZHHyMsFlTOzvg+cD++Y8eidnFxUFCOl3GikN0rkjm+N8veJgQIbuxJ+4GRRLX2bVCz58rkow6YS60c2pjK/tWnyrvMOuk1tLo2hHb9wht80mHNzSX74znkLVwIFgsqnQ7fBx7A94H7L2jbYbEpbD2eQ2p+KQadE23C3NmTs5L3dswgx2rv/35dqYkelu68eaY/BWpP3r+jPYPbBtdFaJIk1ZSyAixf3Ig2Yy+5wo0jvT/hml5DiM+J55Utr3Ao5xAAbfzaMCluEu3821X7FKbjJ0ifOpWSHTsAcAoOJmDCBDwG39Bgq2IA8tKN7F2ZQsK2dBSr/WPZK9BA656hNO8a3CBGTJXJRy3KSS3m0LozJGxLx1JmA8DgoaNtnzBa9Qht0A1Jwd6uI2f+fHI/n4dSYk8cXHteR9Dzz6OLuHBgr13JuTy+cC+n80oB0LgcxxD0MzinAxBpsfC4OoSVtvEsTtKi1aiYPbIDg1oH1V5QkiQ5Tmk+aR9eT3BxPGac0Nz0IZrYO7ApNhYlLuL9Pe9jtNiHF7ih0Q080eEJgt2q98VDCEHR77+T8dZbWFPt1Tou7drhP+EpXOPiajyk2mTMN7F/zSkOrjuD+exnipNeQ0yXINr0DMU31K2OS3hxMvlwMJtV4fjeLA6uO1PeXRbso9nF9o8gJi4IjbZhD6lrKzaSv2ghOfO+wJZjr5N0btWKgAlP4dqt8u5zm5OyuWfeDsw2BW+vXFyDVlCg2guAu03hgfwiLJabmFU8kFIL6DRqPhrVgX4tZVWLJF1Nigrz2fHObfRhu31Bz0nQczKo1WSVZPH+nvdZdmwZAoFeo2dMqzHc1/q+i/aKEUKw91Q+CelFRPoY6NrYF5VKhVJWRu78+WR/+hni7JcjQ9dr8B8/HkP79rUVrkOYy6wc2ZbO/rVnyEszli8PaepFqx4hRMf617uhGWTy4SD5GSUkbEnj8OY0SgrNgL2lcqN2frTuGUpYjHeDvu0H9umu8xZ8Q+6CBSgF9knttJERBDzxBO4DB6JSV55Unc4r4Yb3NlBkyySy8UZyVVsQCNRCcGtRMbcXG5hQ+BAHhH3K+3ZhnrxxcxtahTTMhlWSJF3agi0nKF7+AuOcfrUvaNIPhn8Krr4AxOfEM3PHTHZm7ATAW+/N2NZjub357bg4/d2Go6DEwhOL9rAmMat8WfsILz4Z3ZEAD3uVryUjk+yPPyL/+x/AYh9V1PW6HviNG4ehQ4faCNdhhBCkHsnnwNrTHN+XjVDsH9k6FyeadgqgeddgAht51IvPHpl81CBTiYWjOzNJ2JJGxonC8uUGTx0tu4fQqnsIbt4Nuz0HgDk5mbxvvyVvydLybxC6qCh8H3wQz6FDUGkvXn2kKIIRc1dwwPgDOu9toLLfKuxvLOHRvHyiW90B18/geCEcySgmzNuFViH148UiSZJjKIrg1jmbaXzmJ17XzUcnzOAeAje+D037A/YP1tUpq5m1exbJhckA+Dr7cn+b+7m12a0gtNzx6Vb2nspHp1HTMdKb/afzMZptRPkaWPZIdzwNf783mU+fIXvOxxT8+BPY7O9DLrGx+Iy9F/e+fRvcUO3/VJxXxqGNqSRuSacot6x8uXeQgeZdg2naObBO2xfK5ONfsphtpBzK4djOTE7sy8ZmtQ+Nq1KriGjpQ/OuwTSK9UPTwGcrFIqCccMGcr/5BuP6DeXL9c2b4zfuf7j373/ZF2t2aTZPr5jNztxfUantd4O6lJl5PCeXNmoDDJkFrW92aBySJNVPCemFDHl/I41FMt/7fYpb0XH7ihY3Qt8p4NcUAKti5dfjvzJn3xzOFJ8BwMfZhyBVH7bva46HzptvH+hCqxBPknOM3PnZNs7kl9KvRSCf3d3xgi8y5uRkcubOpeCnZYizd0J0kZF4j7kbzxtvRONWf9tNVIVQBGeO5JGwJZ2k3ZlYzw7fDhAU7UGTjoE07uBf61+MZfJxBSwmG8kHczi2K5Pkg9lYzX9fTN9QV2KuCaZZXCCunvpaK5OjWFJTKfj5Z/J/+BFLSkr5ctfreuAzejSuPXpc9q5EWnEaXxz6gu+PfI9ZsScdMVY1E7LS6FpmgsZ9YNiH4BFyyeNIknR1m/FHAh+vTaKRB/wZuwndjjkgbKBSQ+tboPtTENgSAItiYdmxZXy2/zNSjakACMWJa4MGMqnrg0R72attD6UWMPzDzZhtCnNGd2BQ68obrFqzssj95hvyvltYXo2sMhjwHDwYr9tvx6V1q1r4CziWudTKsd2ZJG5NJ/VYfnl3XbB32W3cMYBG7fzw8HV8d2SZfFSRscBEyqEcTh7IIeVgToXs0d3HmcYd/GkWF4RfuFuDryJQjEYKV66k4Kdl9tlmz152tYcHXjffjPfIO9BFRl72OCcLTvLFoS/4OelnrIoVgKBSZ54tOE3v0hJUOnfoNxU63y+HRZckiTKLjQGz1pOSW8Lw9qG8c50a1dppkPjb3xuFd4FWw6FRTwhogVXYuO3rOSSWLkfjcrp8s06Bnbi56c30i+zHR6tTmL36GEEezvw1oSdu+ot3RVWMRvJ/+JG8hQsxJyWVL3du1QrP4cPxuH4QTr6+Dom/NhnzTSTtyeTYrkzSjhVUWOcT4kpUWz+i2vgR2MjDIZOXyuTjIhRFkHmykOSDOSQfzCErpajCeg8/Zxp3CKBxhwACIt2vioSjeP16ClesoHjd+vK2HACGuDg8b7oJj0EDURsuPa+AEIItaVv45vA3rD+9vnx5a5UPD55JopepCBXYv8UMfAPcZbdZSZL+tvV4DqPmbsOmCCYOiuHhXk0gbR+sfwsSltvvhJzj7EmBR3N+TPXkGKH0GNSc5SV7WJe+FUXYvyC6a90ZEDWIP7cFk5YRxON9Y3iyf7PLlkMIQenOneQtXETRihXlVTJoNLh27YrH4MG49+9XoVrmWGYRaxKyyC81E+XryoCWQRXamdRXxXllJO3OImlPJulJBZz/Se/spiWqtS+97mpeo80HZPJRibSkAn77eD9lxZYKywMi3Ylo7Uujtn74RzT8hMNWUEDx2rUUrliJceNGhMlUvk4bGYHnsGF43jisSlNUGy1GfjvxG98e/pZj+cfKl/d0CeW+lMO0L84HINMthoDhM6Bx7xqPR5Kkq8PXW07y4jL7IGPvjGjHzR3ODo1elA4HlkDSakjeAtbSSvdP9wpjWVAjfhQFnDHnly9XrO5gbMu7Q+6id1QXnNRVG4zLmptL4S+/ULD8N8r27/97hVaLa+dOaLtfxxxzMF+etFbYz13vxOQbmnNnXESD+bwoM1rsX7oPZJN8KBdzqRXfMDfueKFmx0SRyUclyootfP7MBnTOToS38CGqjS8RrXwxeOhq7Bx1QdhslB06RPGGDRg3bKR0/35Q/q4+0kZE4DFwAO4DBuDcuvVlXyxCCHZn7ubHoz+yInkFpWffCAxOBm4yRHDnsZ1EGnMBOKyE851hFM9PeAa9tv6PvidJUt167dd45m48gUat4uv74ujW2K/iBlYzf65dy4rVq2ijO8Od0WXo8o5CXjLnGjMowHZnZ34JjGSNk42is23OANy0bnQO6kzXkK50De5KpEdklRIEc3IyBcuXU7j8twrVMgCn3PzJbNaWkpaxLFcHs+9sp8cRncKYdnNbNA1o+HOwTwGSfqwAq0UhsnXNVjXJ5OMislKK8Al1bdC9VISiYDp2jNJduyjZsQPj5i3Y8vMrbKNv2gT3/v1xHzAAfUxMlRKOw7mH+Sv5L/48+ScpRX83Qo1yj+A2fTDD41fjbrQPNlbm1YSns25guS2Or+/rSvemfhc7tCRJUjlFETyxaC8/70vFx1XHL491J9Tr74aQxSYrvd9aS1aRiRcGt+D+HvYGppiNkLIVjq+FpDWQcQAAC7DFxZmlBj82u2oxaWwVzufj7EOsfyyxAfZHS9+W6DWX7jSQfvAIn838miZJe2mdcxwnoVRYXxwaxWp9GPE+kURd25kXH+yPUwP+TKlJMvm4iihlZZgSEijZvYeSXbso3bkTW0HFhkRqNzdcu3XDtUd33Lp3Rxt8+aGKbYqNA9kH+D7hdzalrSGrLK18ncHJwMDwXgw3CWL3/ojKmGlf4dMYa4+J3Lg2kPiMEm5sF8L7Ixv2KIKSJNWuMouNWz7ezKHUQtqEerJkXFectfYu/dN+P8wn644T5WtgxZM90Tld5EO94DQk/g4Jv8LJjaBYsQEJOh1bfILY4unPHkseFlGxysRJ7URLn5bE+MTQzLtZ+cNNZ2/jkV9i5o5Pt5KQXkSolwuLR7fG/dAeSrZupmT7DkxJJy6Mx+COb8dYXNq2xblNa5ybN8cpMLDBVMnUJJl8NFBKaSllCQmUHYqn7NAhyg4dwpSUVD5YzjkqFxdcYtth6NQJ1y5dcGnX7pKDgJ2TU5rD5tTNrE5ex4bTmzGJvxvcCkULJTH08WjKcy45BCX+CGfnX8Ar0j48ctvbmb32BG+vPIK3QctfT/XE163hdz2WJKl2nc4rYejsjeSVWLiumT/v3h7L7uQ8Hvx6J4qAz+7uRP+qTrtQmk/arl/Y9ecCeqv34Kqyt3MzqeCwfzR7Q1uzV69jb0ESOWU5lR4i1C2UCLdIkpPNuOSX0VpVxig/K43zjqMrPMO5Kh9rmZqSHAMl+d4UZTlhzjSjVi78CFV7eqJv2gTnZjHomzU7+2ja4McXuRyZfNRztmIj5uNJmI4llT+bjh/HcuoUVHI5NL6+uLRpg6FzJwydOuHcsmWVk429mXvZnbmbnRk7ic+Jr7Be2JwRxuY0tYVybV42Q8Q2Wqj/rnJRAlqhvna8vReLRktCeiE3zt6E2abw7u2x3NT+8o1WJUmSKrPteA73fLGDUkvFL1e3dAjj7RHVn+120tL9LNt5jAeCjvJU8EFUR1eA9e9RQIVbIGeCWnDAzZMjKhtHzHkkmrLJsFXewBVAJQTBVhshVivBVhvBVmv5zyFWK0FmKyJPS2mulrzCACg2YM4sBJtS6fGc/P3RRUaijYpEF3nuEYUuIhy1i+PH4XA0mXzUMSEE1qwsLKdPYzl1CvOps89nTmNJOYU1M/Oi+2r8/XBp2QrnVq1wbtUS51atqnQLr8BUQGJuIgm5CSTmJbI/az8nC09esJ2uzJdgozc9hJ7RHjZC8vf+Xa0CWNDyp60jC229yfS7hrdHtKdNmCfpBWXc/ukWknNK6B3jz7x7Ov8nbytKklRz9qTk8fyPB4lPK0SnUTOicxgvDG5ZXg1THekFZfR6aw1lFoVXb2rNXe197FUzB7+HY6tAsVS6X4FazRGdlmStEyecnDnl5kWqXk+KYqJUVL7P+XxsEGoxEWy1EmyzEWKyElIWgK8qBldTEJq0YsxHjl7yfR9A4+eHNigIbXAQTkHBFX8ODkLj64taV787SNSL5OPDDz9k5syZpKen065dO2bPnk1cFaY6rs/Jh7BasRUUYMvPx5qdgzUzA2tmJtbMTCyZmVgzs+y/Z2QgzOZLHkvj74e+cRP00dHomjRGH90YfeNonPz9L7qPVbGSZkwjpTCF5MJkThWdIrkwmaT8pPLRAP+piZMHrU2CFrlZDCjNxq+yjFzrCtE9oWl/RMub+PO4iRd+Okh2sRmNWkX3Jn7sP51PXomFMG8Xfn60Oz6u9ftFIElSwyCEIK/EgrNWjUH373rNfbo+iTd+S8BJreLLsXFc2+RsY3izETIO2ccWyU+BEnv1S6pRsDRJzaEyf9J0UTw3+gauaRJQXq6cshxOFZ0itTiVNGNa+XNacRqpxtTy3oCXohbg6+RChDqQJmU+hOc7EZSr4JFVgiEtH21qNqriksseB0Dt7o6Tjw8aHx80vj44+fii8fG2P3t6oHZ3R+Pu/vezhwdqV9eLTgha0+o8+Vi0aBF33303c+bMoUuXLrz77rssWbKExMREAgICLrmvo5IPoSgohYUoRiM2oxHFaEQpKbE/G889//2wFRZiy8+v8FCKii5/onM0GrTBwWjDwtCFh6ENDUMbHoYuLAxdVBQaT09sig2j1UihqZBC89nH2Z/zTflklmSWPzJKMsgpzcEmbBc9ZahKT3OLlZjiPFqWlRFrMuOpVEw2LM5+aANj7HMq+DaFoDYQcQ04VWy7kWs0M2XZQX7d/3dD1GaBbnw+pjPhPpcelEySJKkuCCEYv3Avv+xLReek5rnrmzP6msgLeqMYTVbeX32UT9cfRwhoHuTOJ3d1JNLXtVrnKjAVkGpMJa04jZVHD/NrfDwqVTbeunS0ukJyVVasVbhD7FYi8C8E30JBULETwUYdAcVqfAoUPAssuOVX3rakSlQq1G5uqN3d0Li5o3Z1Re3qijY0lOCXX7qyY15EnScfXbp0oXPnznzwwQcAKIpCeHg4jz32GJMnT66wrclkwnTeQFiFhYWEh4fXePJxZO1ybOOerpFjlelVlBjUFLvZH0VuaordVBS7qilyU5HnriLXDUxOYEXBohJYEFhUChYEVpXAgoLtCmottEJFoE1DhMVKtLmEaEspURYrMWYzHuf9c5ao3TgiIthvCSVRhOMVGcv/br0eD59LJ3//dOB0AdtO5BDuY6BP8wC0skuZJEn1WJnFxvjv9rAiPgOACB8D3Rr7otWosSoKJWYbG49mk2O0350eGRfB1KFXVtXzTyeyjTz8zW4Op9kHA+kUCIM9dhJUvA6dKYFsjYoMJw0ZGg0FGg15Gi25Gi2FahXFKgXlIp8JKiEwlIFnCXiUgKdR4HHez54lYCgDg0ngagJXk/137cW/q5Lvq6Prpn3/Oubz1WnyYTabMRgMLF26lJtuuql8+ZgxY8jPz2fZsmUVtn/ppZd4+eWXLzhOTScfW1d+g+djr2HWQJkOSvVQqoMyLZTpVJTq7cvLtPZ1RmcVRS5Q5ALFLn//bHQGpYYHlXFWFDzOPtwVBQ9F4GmzEWCzEWC1Pweeffaz2Tj/498iNCSJEBJEOAlKBAkinEQlgjR8ABX+7nomDWrOLR1CZRsNSZL+E2yK4LvtKcz8M5GC0srbbUT6GnjuhhYMbFWz00GUWWx8uOYYn204Ttl584V5UMwAzS56qffRSnWCRuqMCvsJoFilokCjpkCtplitpujs4++fVeU/F5/9/fz1/7zLorUKDKZziYk9OXE2g4sZ3FUwfeahGo29TpOP1NRUQkND2bx5M127di1fPnHiRNatW8e2bdsqbF9bdz6On4rnkz8mgkZln4eEis+VLa2wXFVx+cX2V6vUaFGjRYMTGrSq835Gg9PZdeceOjRoVU7nTcJ23j9P+aLz16mwapwp1gdT4hJIoUsYirpi+4tzZQ31cqZXTECNZPSSJEkNTYnZyqrDmSRlFaMI0KpVOGnUxAS50b2J/8XHEakBOcUmViVkkl5Qhsl64S0Irc2IZ1kaHuZ0PEwZ6GxGNIoZjbBSYWraKhIIzCiUYrU/VGefsWFVKVixPyxnnw06Lybf9WUNRPq36iQfdT4mtl6vR693/FgR0eEtmfHArw4/jyRJklQ/GHRODG0XUifn9nXTM6JTeJ2cuyGo8bTPz88PjUZDRkbFW0oZGRkEBcnZTiVJkiTpv67Gkw+dTkfHjh1ZtWpV+TJFUVi1alWFahhJkiRJkv6bHFLt8tRTTzFmzBg6depEXFwc7777LkajkXvvvdcRp5MkSZIkqQFxSPJx++23k5WVxZQpU0hPTyc2NpY//viDwMAqjtUvSZIkSdJVSw6vLkmSJEnSv1adz285YpQkSZIkSbVKJh+SJEmSJNUqmXxIkiRJklSrZPIhSZIkSVKtksmHJEmSJEm1SiYfkiRJkiTVKpl8SJIkSZJUq2TyIUmSJElSrarzWW3/6dyYZ4WFhXVcEkmSJEmSqurc53ZVxi6td8lHUVERAOHhcipiSZIkSWpoioqK8PT0vOQ29W54dUVRSE1Nxd3dHZVKVaPHLiwsJDw8nFOnTl2VQ7df7fHB1R+jjK/hu9pjlPE1fI6KUQhBUVERISEhqNWXbtVR7+58qNVqwsLCHHoODw+Pq/afCq7++ODqj1HG1/Bd7THK+Bo+R8R4uTse58gGp5IkSZIk1SqZfEiSJEmSVKv+U8mHXq9n6tSp6PX6ui6KQ1zt8cHVH6OMr+G72mOU8TV89SHGetfgVJIkSZKkq9t/6s6HJEmSJEl1TyYfkiRJkiTVKpl8SJIkSZJUq2TyIUmSJElSrZLJhyRJkiRJtapBJx+5ubmMGjUKDw8PvLy8uO+++yguLr7kPp9++im9evXCw8MDlUpFfn7+FR13//799OjRA2dnZ8LDw3nzzTdrMrQql+OfysrKeOSRR/D19cXNzY1bbrmFjIyM8vXz589HpVJV+sjMzARg7dq1la5PT09vEDEClZZ/4cKFFbZZu3YtHTp0QK/X06RJE+bPn1/T4Tkkvn379jFy5EjCw8NxcXGhRYsWvPfeexfE5qhr+OGHHxIVFYWzszNdunRh+/btl9x+yZIlNG/eHGdnZ9q0acNvv/1WYb0QgilTphAcHIyLiwv9+vXj6NGjFba5kr/jlarJ+CwWC5MmTaJNmza4uroSEhLC3XffTWpqaoVjREVFXXCtpk+fXu/jA7jnnnsuKPugQYMqbFOb1w9qPsaLvWfOnDmzfJv6eg0PHTrELbfcUl6+d99994qOWZX33WoRDdigQYNEu3btxNatW8WGDRtEkyZNxMiRIy+5z6xZs8S0adPEtGnTBCDy8vKqfdyCggIRGBgoRo0aJQ4ePCi+++474eLiIj755JM6j2/cuHEiPDxcrFq1SuzcuVNcc801olu3buXrS0pKRFpaWoXHwIEDRc+ePcu3WbNmjQBEYmJihe1sNluNxueoGIUQAhBffPFFhfKXlpaWrz9+/LgwGAziqaeeEvHx8WL27NlCo9GIP/74o97H9/nnn4vx48eLtWvXiqSkJPH1118LFxcXMXv27PJtHHUNFy5cKHQ6nZg3b544dOiQeOCBB4SXl5fIyMiodPtNmzYJjUYj3nzzTREfHy9eeOEFodVqxYEDB8q3mT59uvD09BQ//fST2Ldvn7jxxhtFo0aNKlyvK/k71of48vPzRb9+/cSiRYtEQkKC2LJli4iLixMdO3ascJzIyEjxyiuvVLhWxcXF9T4+IYQYM2aMGDRoUIWy5+bmVjhObV0/R8X4z/fMefPmCZVKJZKSksq3qa/XcPv27eLpp58W3333nQgKChKzZs26omNW5X23Ohps8hEfHy8AsWPHjvJlv//+u1CpVOLMmTOX3f/cm/M/k4+qHPejjz4S3t7ewmQylW8zadIkERMT8y+jql45/ik/P19otVqxZMmS8mWHDx8WgNiyZUul+2RmZgqtViu++uqr8mUX+9vUNEfGCIgff/zxoueeOHGiaNWqVYVlt99+uxg4cOAVRnOh2rqGQgjx8MMPi969e5f/7qhrGBcXJx555JHy3202mwgJCRHTpk2rdPsRI0aIwYMHV1jWpUsX8b///U8IIYSiKCIoKEjMnDmzfH1+fr7Q6/Xiu+++E0L8+9d6ddR0fJXZvn27AERycnL5ssjIyEo/FGqaI+IbM2aMGDZs2EXPWZvXT4jauYbDhg0Tffr0qbCsvl7D812sjJc75pW+L11Kg6122bJlC15eXnTq1Kl8Wb9+/VCr1Wzbts2hx92yZQvXXXcdOp2ufJuBAweSmJhIXl7eFZ+7uuX4p127dmGxWOjXr1/5subNmxMREcGWLVsq3eerr77CYDBw6623XrAuNjaW4OBg+vfvz6ZNm/5lRBdydIyPPPIIfn5+xMXFMW/ePMR54+lt2bKlwjHAfg0v9ne6ErV1DQEKCgrw8fG5YHlNXkOz2cyuXbsqlE2tVtOvX7+Llu1yf+cTJ06Qnp5eYRtPT0+6dOlSvo2jXuu1EV9lCgoKUKlUeHl5VVg+ffp0fH19ad++PTNnzsRqtV55MJVwZHxr164lICCAmJgYHnroIXJyciocozauH9TONczIyGD58uXcd999F6yrj9ewJo55pe9Ll1LvZrWtqvT0dAICAiosc3JywsfH51/Va1fluOnp6TRq1KjCNoGBgeXrvL29r/j81SlHZfvodLoL3tQCAwMvus/nn3/OnXfeiYuLS/my4OBg5syZQ6dOnTCZTMydO5devXqxbds2OnTo8O8C+0d5HRXjK6+8Qp8+fTAYDKxYsYKHH36Y4uJixo8fX36cc9fs/GMUFhZSWlpa4e9RH+M73+bNm1m0aBHLly8vX+aIa5idnY3NZqv075aQkHDReCrb/vzX0rlll9rGEa/1f3JEfP9UVlbGpEmTGDlyZIXZRMePH0+HDh3w8fFh8+bNPPvss6SlpfHOO+/8y6j+5qj4Bg0axM0330yjRo1ISkriueee4/rrr2fLli1oNJpau35QO9fwyy+/xN3dnZtvvrnC8vp6DWvimFfyvnQ59S75mDx5MjNmzLjkNocPH66l0tS8+hTfli1bOHz4MF9//XWF5TExMcTExJT/3q1bN5KSkpg1a9YF21amPsT44osvlv/cvn17jEYjM2fOLE8+/o36EN85Bw8eZNiwYUydOpUBAwaUL/+311CqeRaLhREjRiCE4OOPP66w7qmnnir/uW3btuh0Ov73v/8xbdq0ej/HyB133FH+c5s2bWjbti2NGzdm7dq19O3btw5L5hjz5s1j1KhRODs7V1jekK9hXah3yceECRO45557LrlNdHQ0QUFB5b0zzrFareTm5hIUFHTF56/KcYOCgi5o5Xvu98ud25HxBQUFYTabyc/Pr5ChZmRkVLrP3LlziY2NpWPHjpcsD0BcXBwbN2687HZQv2I8p0uXLrz66quYTCb0ev1Fr6GHh8dl73rUl/ji4+Pp27cvDz74IC+88MIlywPVu4aV8fPzQ6PRVPp3u1Q8l9r+3HNGRgbBwcEVtomNjS3fxhGv9X9yRHznnEs8kpOTWb16dYW7HpXp0qULVquVkydPVkgi/w1Hxne+6Oho/Pz8OHbsGH379q216weOj3HDhg0kJiayaNGiy5alvlzDmjjmlb7vXtIVtRSpB841Ytq5c2f5sj///LPGGpxe6rjnGpyazebybZ599lmHNDitTnznGgUtXbq0fFlCQkKljYKKioqEm5tbhR4Sl9KvXz8xfPjwK4jk4hwd4/lee+014e3tXf77xIkTRevWrStsM3LkSIc0OHVEfAcPHhQBAQHimWeeqXJ5auIaxsXFiUcffbT8d5vNJkJDQy/ZmG/IkCEVlnXt2vWCBqdvvfVW+fqCgoJKG5xe6Wu9Omo6PiGEMJvN4qabbhKtWrUSmZmZVSrHggULhFqtvqDXyL/liPj+6dSpU0KlUolly5YJIWr3+gnh2BjHjBlzQU+li6kv1/B8l2pweqljXun77qU02ORDCHv3rfbt24tt27aJjRs3iqZNm1bovnX69GkRExMjtm3bVr4sLS1N7NmzR3z22WcCEOvXrxd79uwROTk5VT5ufn6+CAwMFHfddZc4ePCgWLhwoTAYDA7palvd+MaNGyciIiLE6tWrxc6dO0XXrl1F165dLzj23LlzhbOzc6W9IWbNmiV++ukncfToUXHgwAHx+OOPC7VaLf76668ajc9RMf7888/is88+EwcOHBBHjx4VH330kTAYDGLKlCnl25zravvMM8+Iw4cPiw8//NBhXW1rOr4DBw4If39/MXr06Ard+s7/YHPUNVy4cKHQ6/Vi/vz5Ij4+Xjz44IPCy8tLpKenCyGEuOuuu8TkyZPLt9+0aZNwcnISb731ljh8+LCYOnVqpV1tvby8xLJly8T+/fvFsGHDKu1qe6m/Y02p6fjMZrO48cYbRVhYmNi7d2+F63Wut9zmzZvFrFmzxN69e0VSUpJYsGCB8Pf3F3fffXe9j6+oqEg8/fTTYsuWLeLEiRPir7/+Eh06dBBNmzYVZWVl5ceprevniBjPKSgoEAaDQXz88ccXnLM+X0OTyST27Nkj9uzZI4KDg8XTTz8t9uzZI44ePVrlYwpR9c+WqmrQyUdOTo4YOXKkcHNzEx4eHuLee+8VRUVF5etPnDghALFmzZryZVOnThXABY8vvviiyscVQoh9+/aJ7t27C71eL0JDQ8X06dPrRXylpaXi4YcfFt7e3sJgMIjhw4eLtLS0C47dtWtXceedd1Z63hkzZojGjRsLZ2dn4ePjI3r16iVWr15d4/EJ4ZgYf//9dxEbGyvc3NyEq6uraNeunZgzZ84FY1ysWbNGxMbGCp1OJ6Kjoyv8D9Tn+C72PxwZGVm+jSOv4ezZs0VERITQ6XQiLi5ObN26tXxdz549xZgxYypsv3jxYtGsWTOh0+lEq1atxPLlyyusVxRFvPjiiyIwMFDo9XrRt29fkZiYWGGbqrwma0pNxnfu+lb2OHfNd+3aJbp06SI8PT2Fs7OzaNGihXjjjTcqfHjX1/hKSkrEgAEDhL+/v9BqtSIyMlI88MADFT60hKjd61fTMZ7zySefCBcXF5Gfn3/Buvp8DS/2P3j+2E6XO6YQVf9sqSqVEOf1P5QkSZIkSXKwBjvOhyRJkiRJDZNMPiRJkiRJqlUy+ZAkSZIkqVbJ5EOSJEmSpFolkw9JkiRJkmqVTD4kSZIkSapVMvmQJEmSJKlWyeRDkiRJkqRaJZMPSZIkSZJqlUw+JEmSJEmqVTL5kCRJkiSpVv0f0Q1bBZcWbfUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we test several values of h to see which is the most suitable\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "sd = round(calculate_std_dev(returns_2015_2016), 4)\n",
    "\n",
    "h = [0.005, 0.01, sd, 0.05, 0.1]\n",
    "x = np.linspace(-0.1, 0.1, 1000)\n",
    "for h_i in h:\n",
    "    y = [kernel_density_estimate(xi, returns_2015_2016, h_i) for xi in x]\n",
    "    plt.plot(x, y, label=f\"h={h_i}\")\n",
    "plt.legend()\n",
    "plt.title(\"Kernel density estimate\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated VaR for 2015-2016 :\n",
      "Confiance\tVaR\n",
      "90.0%\t\t-0.027826494540893393\n",
      "95.0%\t\t-0.038123774394452845\n",
      "99.0%\t\t-0.05683570780710162\n"
     ]
    }
   ],
   "source": [
    "def estimate_var(data, h, probability):\n",
    "    # Estimation of the distribution function using the trapezoid method\n",
    "    def cumulative_distribution(x, data, h):\n",
    "        return trapezes(kernel_density_estimate, min(data), x, 1000, data, h)\n",
    "\n",
    "    # Find the value where the distribution function reaches the probability level\n",
    "    lower_bound = min(data)\n",
    "    upper_bound = max(data)\n",
    "\n",
    "    # Dichotomy method to find the value corresponding to the probability level\n",
    "    while upper_bound - lower_bound > 1e-5:\n",
    "        middle = (lower_bound + upper_bound) / 2\n",
    "        if cumulative_distribution(middle, data, h) < probability:\n",
    "            lower_bound = middle\n",
    "        else:\n",
    "            upper_bound = middle\n",
    "\n",
    "    return upper_bound\n",
    "\n",
    "confidence_level = [0.1, 0.05, 0.01]\n",
    "h = 0.01 \n",
    "print(\"Estimated VaR for 2015-2016 :\") \n",
    "print(\"Confiance\\tVaR\")\n",
    "for confidence in confidence_level:\n",
    "    var_2015_2016_kernel = estimate_var(returns_2015_2016, h, confidence)\n",
    "    print(f\"{(1-confidence) * 100}%\\t\\t{var_2015_2016_kernel}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The proportion of returns exceeding VaR is : 1.5717092337917484%\n",
      "There are 8 on 509\n"
     ]
    }
   ],
   "source": [
    "var_2015_2016_kernel = estimate_var(returns_2015_2016, h, 0.05)\n",
    "\n",
    "# Calculating returns between 2017-2018 \n",
    "returns_2017_2018 = calculate_returns(Prix_2017_2018)\n",
    "\n",
    "# proportion of returns exceeding calculated VaR\n",
    "rendements_sup_var = len([x for x in returns_2017_2018 if x <= var_2015_2016_kernel])\n",
    "proportion = rendements_sup_var / len(returns_2017_2018)\n",
    "\n",
    "# Display of the proportion of returns exceeding VaR\n",
    "print(f\"The proportion of returns exceeding VaR is : {proportion*100}%\")\n",
    "print(f\"There are {rendements_sup_var} on {len(returns_2017_2018)}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question B (Ex2, Q5 of TD2) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Expected Shortfall for VaR is : -0.05336181815045533\n"
     ]
    }
   ],
   "source": [
    "# Identify returns less than VaR during 2015-2016\n",
    "losses = [x for x in returns_2015_2016 if x <= var_2015_2016_kernel]\n",
    "\n",
    "# Calculate Expected Shortfall (ES = average loss above VaR)\n",
    "expected_shortfall = sum(losses) / len(losses) \n",
    "\n",
    "print(f\"The Expected Shortfall for VaR is : {expected_shortfall}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question C (Ex2, Q1 and Q2 of TD3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pikands Index :\n",
      "Losses\t\t\tGain\n",
      "0.5089715779341741\t0.5772338569463286\n"
     ]
    }
   ],
   "source": [
    "# retrieves all available prices\n",
    "Prix = []\n",
    "for year in range(2015, 2019):\n",
    "    for month in range(1, 13):\n",
    "        Prix += price_data[year][month]\n",
    "\n",
    "# Calculation of corresponding returns\n",
    "returns = calculate_returns(Prix)\n",
    "\n",
    "sorted_returns = sorted(returns)\n",
    "\n",
    "\n",
    "def k(n):\n",
    "    return np.log(n)\n",
    "\n",
    "def Pickands_estim(data):\n",
    "    n = len(data)\n",
    "    etam = data[n-int(k(n))]-data[n-2*int(k(n))]\n",
    "    den = data[n-2*int(k(n))]-data[n-4*int(k(n))]\n",
    "    res = np.log(etam/den)\n",
    "    return res/np.log(2)\n",
    "\n",
    "\n",
    "# invert everything to be able to calculate the Pikands index on losses\n",
    "sorted_returns_losses = [-x for x in sorted_returns]\n",
    "sorted_returns_losses = sorted_returns_losses[::-1]\n",
    "\n",
    "# Calculating the Pikands index\n",
    "pikands_gain = Pickands_estim(sorted_returns)\n",
    "pikands_losses = Pickands_estim(sorted_returns_losses)\n",
    "\n",
    "print(\"Pikands Index :\")\n",
    "print(\"Losses\\t\\t\\tGain\")\n",
    "print(f\"{-pikands_losses}\\t{pikands_gain}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VaR EVT :\n",
      "Confiance\tVaR loss\t\tVaR gain\n",
      "90.0%\t\t-0.018676315001953256\t0.034136255033759964\n",
      "95.0%\t\t-0.03459467977312686\t0.03787128303972365\n",
      "99.0%\t\t-0.05563667261604414\t0.05522424941671052\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def VaR_EVT(p, data, pickands):\n",
    "    n = len(data)\n",
    "    num = ((k(n)/(n*(1-p)))**pickands - 1) * (data[n-int(k(n))]-data[n-2*int(k(n))])\n",
    "    den = 1 - 2**(-pickands)\n",
    "    return num/den + data[n-int(k(n))]\n",
    "\n",
    "alpha = [0.1, 0.05, 0.01]\n",
    "print(\"VaR EVT :\")\n",
    "print(\"Confiance\\tVaR loss\\t\\tVaR gain\")\n",
    "for p in alpha:\n",
    "    loss = VaR_EVT(1-p, sorted_returns_losses, pikands_losses)\n",
    "    gain = VaR_EVT(1-p, sorted_returns, pikands_gain)\n",
    "    print(f\"{(1-p)*100}%\\t\\t{-loss}\\t{gain}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question D (Ex2, Q3 and Q4 of TD4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import pandas as pd\n",
    "import numpy.random as rd\n",
    "from scipy.stats import linregress\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "data=pd.read_excel(\"Dataset TD4.xlsx\")\n",
    "nouveaux_noms = {\n",
    "    'transaction date (1=1day=24 hours)': 'transaction_date',\n",
    "    'bid-ask spread': 'bid-ask_spread',\n",
    "    'volume of the transaction (if known)': 'volume',\n",
    "    'Sign of the transaction': 'sign',\n",
    "    'Price (before transaction)': 'price'\n",
    "}\n",
    "\n",
    "data = data.rename(columns=nouveaux_noms) #We rename the columns to use them more easily\n",
    "volume = []\n",
    "price = []\n",
    "sign = []\n",
    "for i in range(len(data) - 1):\n",
    "    if pd.notna(data['volume'].iloc[i]):\n",
    "        volume.append(data['volume'].iloc[i])\n",
    "        price.extend([data['price'].iloc[i], data['price'].iloc[i + 1]]) #We retrieve the price at which the transaction was carried out and the price at t+1 to see the impact of the order on the market.\n",
    "        sign.append(data['sign'].iloc[i])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Annualized volatility: 0.05849549575608429\n",
      "\n",
      "Estimation of eta: 2.280706114725874e-08\n",
      "Mean Squared Error: 0.00023873091719548648 ; R^2: 0.9607480862831544\n",
      "\n",
      "Estimation of gamma: 0.00050237805904092\n",
      "Mean Squared Error: 0.0004907505222142958 ; R^2: 0.9193112589657632\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAls0lEQVR4nO3deXxV9Z3/8dcn+0b2hZAESCCALKISFhcsrUXR+hOdsZbWVtqqTK1t7dRO1fY3nd9MpzN27Kat2iLaautGrVZaxQ1FHWUxyBLZQwIkEJJAWAJhS/L9/XEPNECEkO3cm/t+Ph553HO/55zk833A477vOd9zvsecc4iIiET4XYCIiAQHBYKIiAAKBBER8SgQREQEUCCIiIgnyu8COiszM9MNHjzY7zJERELKsmXLdjrnstpbF7KBMHjwYEpLS/0uQ0QkpJjZlo9bp1NGIiICKBBERMSjQBAREUCBICIiHgWCiIgACgQREfEoEEREBAjDQFhVvYefvLIOTfstInKisAuEFVV7eHjhJpZt2e13KSIiQSXsAuH6cfmkxEfzyLsVfpciIhJUwi4QEmKiuHHiQF5bU8vmnQf8LkdEJGiEXSAAzLxoMFERxu/eq/S7FBGRoBGWgZCTHMc1Y/OYW1rNnqYjfpcjIhIUwjIQAG6ZXMjBoy08uWSr36WIiASFsA2Ec3KTmVycyePvb+ZIc6vf5YiI+C5sAwHglslF1DUeZt7K7X6XIiLiu7AOhEuLMxmWk8Scdyt0o5qIhL0zBoKZPWZmdWb2UZu2+8xsnZmtMrMXzCy1zbp7zKzczNab2RVt2seZWZm37gEzM6891sye9dqXmNng7u3iafvGLZcUsW5HI++V7+qtPysiEpQ6coTwe2DaSW2vA6Odc+cCG4B7AMxsJDADGOXt85CZRXr7PAzMAoq9n2O/82Zgt3NuKPAL4Ced7UxnTD9/AJlJsbpRTUTC3hkDwTn3DtBwUttrzrlm7+1iIN9bng4845w77JyrBMqBCWaWCyQ75xa5wLmZJ4Br2+zzuLf8HHDZsaOH3hAbFcnMCwfx9oZ6NtQ29tafFREJOt0xhvBVYL63nAdUtVlX7bXlecsnt5+wjxcye4GM9v6Qmc0ys1IzK62vr++G0gO+OGkQcdERzNFRgoiEsS4Fgpn9AGgGnjzW1M5m7jTtp9vn1EbnZjvnSpxzJVlZWWdb7sdKS4zh+nH5/GX5duoaD3Xb7xURCSWdDgQzmwlcDdzo/n6JTjVQ0GazfGC7157fTvsJ+5hZFJDCSaeoesPNlxRxtLWVPy7a0tt/WkQkKHQqEMxsGnAXcI1zrqnNqnnADO/KoUICg8dLnXM1QKOZTfLGB24CXmyzz0xv+XrgTefDNaCFmYl8+pwc/rB4CwePtPT2nxcR8V1HLjt9GlgEDDezajO7Gfg10A943cxWmNlvAJxzq4G5wBrgFeB259yxT9fbgDkEBpo38fdxh0eBDDMrB74D3N1dnTtbt1xSyO6mo/z5w+ozbywi0sdYqN6QVVJS4kpLS7v1dzrnmP7ge+w/1Mwb3/kEERG9drGTiEivMLNlzrmS9taF9Z3KJzMzbplcRMXOAyxYV+d3OSIivUqBcJKrRvcnLzVeN6qJSNhRIJwkKjKCr1w8mKWVDayq3uN3OSIivUaB0I7PjS+gX2wUc97VE9VEJHwoENrRLy6aGRMKeKmshm17DvpdjohIr1AgfIwvX1wIwO/13GURCRMKhI+RlxrPVWNyeWZpFY2HjvpdjohIj1MgnMatkwtpPNzMsx9UnXljEZEQp0A4jXPzU5lQmM7v3ttMc4ueuywifZsC4QxunVzEtj0Hmf/RDr9LERHpUQqEM7hsRDZFmYl67rKI9HkKhDOIiDC+ekkhK6v3sqhCz10Wkb5LgdAB14/Lp39yHD97bYOOEkSkz1IgdEBcdCTfuqyYZVt286YmvRORPkqB0EGfLclncEYC9726ntZWHSWISN+jQOig6MgIvnP5cNbtaOSvq7afeQcRkRCjQDgLV4/JZWRuMj9/fQNHdV+CiPQxCoSzEBFh/MsVw9myq0l3L4tIn6NAOEtThmcxfnAaDyzYyMEjLWfeQUQkRCgQzpKZ8S9XjKCu8TCPL9rsdzkiIt1GgdAJEwrTmTI8i4cXbmLvQc2EKiJ9gwKhk757+XD2HjzKI+/o2csi0jecMRDM7DEzqzOzj9q0pZvZ62a20XtNa7PuHjMrN7P1ZnZFm/ZxZlbmrXvAzMxrjzWzZ732JWY2uJv72CNG56Vw9bm5PPZeJfWNh/0uR0SkyzpyhPB7YNpJbXcDC5xzxcAC7z1mNhKYAYzy9nnIzCK9fR4GZgHF3s+x33kzsNs5NxT4BfCTznamt915+XAON7fy4FvlfpciItJlZwwE59w7QMNJzdOBx73lx4Fr27Q/45w77JyrBMqBCWaWCyQ75xa5wGRAT5y0z7Hf9Rxw2bGjh2BXmJnIDSX5PLlkC1UNTX6XIyLSJZ0dQ8hxztUAeK/ZXnse0PYC/WqvLc9bPrn9hH2cc83AXiCjvT9qZrPMrNTMSuvr6ztZevf61mXFmBn3L9jodykiIl3S3YPK7X2zd6dpP90+pzY6N9s5V+KcK8nKyupkid0rNyWemRcO4vkPq9lY2+h3OSIindbZQKj1TgPhvR6bArQaKGizXT6w3WvPb6f9hH3MLApI4dRTVEHttilDSYiJ4mevbfC7FBGRTutsIMwDZnrLM4EX27TP8K4cKiQweLzUO63UaGaTvPGBm07a59jvuh5404XYQwfSE2O4dXIRr6zewcqqPX6XIyLSKR257PRpYBEw3Myqzexm4F5gqpltBKZ673HOrQbmAmuAV4DbnXPH5ne4DZhDYKB5EzDfa38UyDCzcuA7eFcshZqbJxeSkRjDfa+u97sUEZFOsRD7Mn5cSUmJKy0t9buMEzz6v5X86G9reOqWiVw0NNPvckRETmFmy5xzJe2t053K3ejGiQMZkBLHT15dr0dtikjIUSB0o7joSO74dDErq/bw2ppav8sRETkrCoRu9o8X5FOUlchPX11Pix61KSIhRIHQzaIiI7hz6nA21u3nL8u3+V2OiEiHKRB6wJWj+zM6L5lfvLGBI8161KaIhAYFQg8IPGpzBNW7D/LHxVv8LkdEpEMUCD3k0uJMJhdn8vPXN1C775Df5YiInJECoYeYGT+aPpojLa38x9/W+F2OiMgZKRB60ODMRL75yaG8tKqGhevrzryDiIiPFAg9bNYniijKSuSHL67m0NGWM+8gIuITBUIPi42K5MfXjmFrQxO/flNPVhOR4KVA6AUXDsngHy7I47fvbKK8Ts9MEJHgpEDoJT+46hwSYqL4/gsfaZ4jEQlKCoRekpEUyz1XjmBpZQPPLas+8w4iIr1MgdCLbigpoGRQGv/18loaDhzxuxwRkRMoEHpRRITx4+vG0HiomXvnr/W7HBGREygQetnw/v24ZXIRc0urWVoZUo+OFpE+ToHgg29dNpS81Hh+8EKZJr8TkaChQPBBQkwUP7p2FBvr9vPIuxV+lyMiAigQfPOpETlMG9WfBxZsZOuuJr/LERFRIPjp364ZSVSE8cN5ujdBRPynQPBRbko837l8OAvX1zP/ox1+lyMiYa5LgWBm/2xmq83sIzN72szizCzdzF43s43ea1qb7e8xs3IzW29mV7RpH2dmZd66B8zMulJXKJl54SBGDUjm3/+6msZDR/0uR0TCWKcDwczygG8BJc650UAkMAO4G1jgnCsGFnjvMbOR3vpRwDTgITOL9H7dw8AsoNj7mdbZukJNVGQEP75uDHWNh/nZaxv8LkdEwlhXTxlFAfFmFgUkANuB6cDj3vrHgWu95enAM865w865SqAcmGBmuUCyc26RC5xIf6LNPmHhvIJUvjRpEE8s2kxZ9V6/yxGRMNXpQHDObQN+CmwFaoC9zrnXgBznXI23TQ2Q7e2SB1S1+RXVXluet3xy+ynMbJaZlZpZaX19fWdLD0rfvWI4GUmxfP+FMlpaNcAsIr2vK6eM0gh86y8EBgCJZvbF0+3STps7Tfupjc7Nds6VOOdKsrKyzrbkoJYcF80Prx5J2ba9PLFos9/liEgY6sopo08Dlc65eufcUeB54CKg1jsNhPd67NmR1UBBm/3zCZxiqvaWT24PO1efm8snhmXxk1fW6bkJItLruhIIW4FJZpbgXRV0GbAWmAfM9LaZCbzoLc8DZphZrJkVEhg8XuqdVmo0s0ne77mpzT5hxcy47/pzSYyJ4htPLdcjN0WkV3VlDGEJ8BzwIVDm/a7ZwL3AVDPbCEz13uOcWw3MBdYArwC3O+eOfeLdBswhMNC8CZjf2bpCXXZyHD/97FjW7Wjkv1/WjKgi0nssVO+QLSkpcaWlpX6X0WN+9Lc1PPq/lTxyUwlTR+b4XY6I9BFmtsw5V9LeOt2pHKS+N204owYk8y/PrWTH3kN+lyMiYUCBEKRioyL51efP50hzK99+drkuRRWRHqdACGJFWUn8+zWjWFzRwMMLy/0uR0T6OAVCkLt+XD7XjB3AL97YyLItesKaiPQcBUKQMzN+fN1oBqTG8a2nV7D3oCbAE5GeoUAIAf3ionlgxvnU7jvE958v07MTRKRHKBBCxPkD07jz8uG8VFbDsx9UnXkHEZGzpEAIIf90aRGXDM3k//11taa2EJFup0AIIRERxs9vGEuCprYQkR6gQAgx2clx/Myb2uLe+ev8LkdE+hAFQgj65IhsvnpxIb9/fzNvrKn1uxwR6SMUCCHqriuHMzJXU1uISPdRIISo2KhIfvWF8zmsqS1EpJsoEELYkDZTWzz0lqa2EJGuUSCEuL9PbbGBN9dpPEFEOk+BEOLMjP/+hzGck5vMN59azprt+/wuSURClAKhD0iMjeLRmePpFxfNzY9/QO0+DTKLyNlTIPQR/VPiePTLJew9eJRbHi+l6Uiz3yWJSIhRIPQhowak8KvPn8/q7Xv59jMrdOWRiJwVBUIfc9k5Ofzr1SN5bU0t985f63c5IhJCovwuQLrfVy4uZPPOAzzybiWDMxO5ceIgv0sSkRCgQOij/vXqkWxpaOKHL66mIC2BS4dl+V2SiAS5Lp0yMrNUM3vOzNaZ2Vozu9DM0s3sdTPb6L2mtdn+HjMrN7P1ZnZFm/ZxZlbmrXvAzKwrdQlERUbw6y9cQHF2Erc/+SEbajVdtoicXlfHEO4HXnHOjQDGAmuBu4EFzrliYIH3HjMbCcwARgHTgIfMLNL7PQ8Ds4Bi72daF+sSICk2ike/PJ64mEi+8rsPqG887HdJIhLEOh0IZpYMXAo8CuCcO+Kc2wNMBx73NnscuNZbng4845w77JyrBMqBCWaWCyQ75xa5wLMhn2izj3RRXmo8j84sYdeBw9z6RKmeoSAiH6srRwhFQD3wOzNbbmZzzCwRyHHO1QB4r9ne9nlA22c/Vntted7yye2nMLNZZlZqZqX19fVdKD28nJufyv0zzmdl9R7unLuSVl2OKiLt6EogRAEXAA87584HDuCdHvoY7Y0LuNO0n9ro3GznXIlzriQrS4OkZ+OKUf2558oRvFRWw09fW+93OSIShLoSCNVAtXNuiff+OQIBUeudBsJ7rWuzfUGb/fOB7V57fjvt0s1unVzE5ycM5KGFm5hbWnXmHUQkrHQ6EJxzO4AqMxvuNV0GrAHmATO9tpnAi97yPGCGmcWaWSGBweOl3mmlRjOb5F1ddFObfaQbmRn/MX0Uk4sz+f7zZby/aaffJYlIEOnqVUbfBJ40s1XAecB/AfcCU81sIzDVe49zbjUwl0BovALc7pw7NsJ5GzCHwEDzJmB+F+uSjxEdGcGDN15AYWYiX/vDMtbWaHZUEQmwwIU9oaekpMSVlpb6XUbIqmpo4obfLuLQ0RaeunUS5+Qm+12SiPQCM1vmnCtpb53mMgpTBekJPH3rJOKiI/nCI4v1HAURUSCEs8GZiTwzKxAKN85RKIiEOwVCmBuUEQiF+OhIvqBQEAlrCgRhUEYiT8+aRIIXCqu37/W7JBHxgQJBgGNHCheSEB3JjXOW8NE2hYJIuFEgyHEDMxKOh8IXH1UoiIQbBYKc4FgoJMZE6UhBJMwoEOQUgVCYRFKsQkEknCgQpF0F6QoFkXCjQJCPpVAQCS8KBDkthYJI+FAgyBm1DYUvPLKY5Vt3+12SiPQABYJ0yLFQSE2IYcbsxby0qsbvkkSkmykQpMMK0hN44esXMTovhduf+pAH3yonVGfLFZFTKRDkrGQkxfLkLROZft4A7nt1Pd/90yqONLf6XZaIdIMovwuQ0BMXHckvP3ceRZlJ/OKNDVTtbuK3XxxHWmKM36WJSBfoCEE6xcy449PF3D/jPFZs3cN1D71HRf1+v8sSkS5QIEiXTD8vj6dunci+Q81c99D7LNq0y++SRKSTFAjSZSWD0/nL1y8mq18sNz22hLmlVX6XJCKdoECQbjEwI4E/33YREwsz+N5zq/jJK+tobdUVSCKhRIEg3SYlPprffWU8n58wkIcXbuL2pz7k4JEWv8sSkQ5SIEi3io6M4L+uG83//cw5vLJ6BzNmL6Ju3yG/yxKRDuhyIJhZpJktN7O/ee/Tzex1M9vovaa12fYeMys3s/VmdkWb9nFmVuate8DMrKt1iX/MjFsmFzH7SyVsqN3PtQ++p2c1i4SA7jhCuANY2+b93cAC51wxsMB7j5mNBGYAo4BpwENmFunt8zAwCyj2fqZ1Q13is6kjc/jT1y6k1cG1D73HE4s2685mkSDWpUAws3zgM8CcNs3Tgce95ceBa9u0P+OcO+ycqwTKgQlmlgskO+cWucCnxRNt9pEQNzovhb996xIuGpLBD19cza1PLGP3gSN+lyUi7ejqEcIvge8BbecuyHHO1QB4r9leex7Q9nrEaq8tz1s+uf0UZjbLzErNrLS+vr6LpUtvyUyK5bGZ4/nXq0fy9oY6pt3/Du9v2ul3WSJykk4HgpldDdQ555Z1dJd22txp2k9tdG62c67EOVeSlZXVwT8rwSAiwrj5kkJe+PrFx5/XfN+r6zjaonmQRIJFV44QLgauMbPNwDPAp8zsj0CtdxoI77XO274aKGizfz6w3WvPb6dd+qDReSn89ZuX8Nlx+Tz41iZu+O0iqhqa/C5LROhCIDjn7nHO5TvnBhMYLH7TOfdFYB4w09tsJvCitzwPmGFmsWZWSGDweKl3WqnRzCZ5Vxfd1GYf6YMSY6P4n+vH8qvPn0957X6uuv9dXlyxze+yRMJeT9yHcC8w1cw2AlO99zjnVgNzgTXAK8Dtzrljdy3dRmBguhzYBMzvgbokyPyfsQN4+Y7JFOckccczK/jun1Zy4HCz32WJhC0L1csAS0pKXGlpqd9lSDdobmnl/gUb+fVb5QzOSORXnz+f0Xkpfpcl0ieZ2TLnXEl763SnsvguKjKCOy8fzlO3TOLgkRaue+g95rxbobmQRHqZAkGCxoVDMph/x2SmDM/mP19ay8zfLWXrLg04i/QWBYIElbTEGGZ/aRz/ee1oPtyym6m/eJsH3yrXYzpFeoECQYKOmfHFSYN4485P8KkR2dz36nqueuBdFlfo4TsiPUmBIEErNyWeh784jse+XMKhoy3MmL2YO+euZNf+w36XJtInKRAk6H1qRA6v//MnuG3KEF5csY3Lfv42z36wVYPOIt1MgSAhIT4mkrumjQjct5CdxF1/LuNzsxexfkej36WJ9BkKBAkpw3L68eysC/mf68+lvG4/n3ngXe6dv46mI7qhTaSrFAgSciIijBtKClhw5xSuOz+P37y9iak/f4c319X6XZpISFMgSMhKT4zhvs+O5dlZk4iPieSrvy/ln/5QSuXOA36XJhKSFAgS8iYWZfDytybzvWnDeWfDTj7987e55/kyduzVs5xFzobmMpI+pa7xEA++Wc5TS7cSYcbMiwZz2yeGkJYY43dpIkHhdHMZKRCkT6pqaOIXb2zgheXbSIqJYtalRXz1kkISY6P8Lk3EVwoECVvrdzTy09fW8/qaWjKTYrj9k0P5wsSBxEZF+l2aiC8UCBL2Pty6m/teWc+iil3kpcbzz1OHcd35eURGtPcEV5G+S9NfS9i7YGAaT906kT/cPIH0xBi++6eVTPvlO7zy0Q5C9UuRSHdTIEjYMDMmF2cx7xsX8/CNF9DqHF/74zKufeh93lhTq6kwJOzplJGEreaWVp5fvo3739jItj0HKcpK5JZLiviHC/KIi9YYg/RNGkMQOY3mllZe/mgHj7xTQdm2vWQkxvClCwfxpUmDyEiK9bs8kW6lQBDpAOccSyobeOSdChasqyM2KoJ/HJfPzZcUMiQrye/yRLrF6QJBF2WLeMyMSUUZTCrKoLxuP4/+bwXPLavm6aVbuWxEDrMuLWL84DTMdGWS9E06QhA5jZ37D/PEoi38YdFmdjcdZWx+CrdeWsS0Uf2JitQ1GRJ6euSyUzMrMLO3zGytma02szu89nQze93MNnqvaW32ucfMys1svZld0aZ9nJmVeeseMH0FkyCRmRTLd6YO4/27L+M/rx3NvkPNfOOp5Uz56UJ++/Ym6hv19DbpOzp9hGBmuUCuc+5DM+sHLAOuBb4MNDjn7jWzu4E059xdZjYSeBqYAAwA3gCGOedazGwpcAewGHgZeMA5N/90f19HCOKH1lbHG2trmfNuJUs3NxAVYXxqRDafG1/AJ4Zl6ahBgl6PjCE452qAGm+50czWAnnAdGCKt9njwELgLq/9GefcYaDSzMqBCWa2GUh2zi3yin2CQLCcNhBE/BARYVw+qj+Xj+pPed1+/lRaxZ8/rOa1NbXkJMfyjxfkc0NJAYMzE/0uVeSsdcugspkNBs4HlgA5XljgnKsxs2xvszwCRwDHVHttR73lk9vb+zuzgFkAAwcO7I7SRTptaHYS91x1Dt+9Yjhvrqtj7gdV/ObtTTy0cBMTC9P53PgCrhydS3yM7mmQ0NDlQDCzJODPwLedc/tOc/q/vRXuNO2nNjo3G5gNgVNGZ1+tSPeLjozgilH9uWJUf2r3HeK5ZdX8qbSK78xdyb+9uJprzhvA58YXMCYvRVcoSVDrUiCYWTSBMHjSOfe811xrZrne0UEuUOe1VwMFbXbPB7Z77fnttIuEnJzkOG7/5FC+PmUISyobmPtB4JTSk0u2MqJ/P64fl89VY3IZkBrvd6kip+jKoLIRGCNocM59u037fcCuNoPK6c6575nZKOAp/j6ovAAo9gaVPwC+SeCU08vAr5xzL5/u72tQWULFvkNH+evK7cz9oIqV1XsBOH9gKp8Zk8uVY3LJUzhIL+qRO5XN7BLgXaAMaPWav0/gQ30uMBDYCnzWOdfg7fMD4KtAM4FTTPO99hLg90A8gcHkb7ozFKZAkFBUufMAL5fV8HJZDau37wPgvIJj4dCf/LQEnyuUvk5TV4gEoc07D/DSSeEwtiCVz4zpz5WjcylIVzhI91MgiAS5Lbv+Hg4fbfPCIT+Fq8bkctUYhYN0HwWCSAjZsusAL5ft4OWyGsq2BcYchmYnMWVYFlOGZzO+ME2PAJVOUyCIhKitu5p4bc0O3t5Qz5KKBo60tJIQE8lFQzKZMjyLKcOzNO4gZ0WBINIHHDjczKJNu1i4oY6F6+up3n0QgOLsJC8cshk/OJ2YKE2fIR9PgSDSxzjn2FS/n4Xr61m4vp6llYGjh8SYSC4aGjh6uGhIJoMzEnQznJxAz0MQ6WPMjKHZ/Ria3Y9bJhcdP3p4a33g6OH1NbUAZPeLZVJRBhOL0plYmMGQrEQFhHwsBYJIH5AYG8WnR+bw6ZE5OOeo2HmAxRW7WFLRwOKKXcxbGbj5PzMplolF6UwqTGdSUQZDs5MUEHKcAkGkjzEzhmQlMSQriRsnDsI5x+ZdTSyp2BUIicoGXlpVA0BGYgwTvHCYWJTOsOx+REQoIMKVAkGkjzMzCjMTKcxMZMaEgTjnqGo4yOKKXSyuDBxFzP9oBwD94qIYm5/KeQXez8BUMpNife6B9BYFgkiYMTMGZiQwMCOBG8YH5pusamhiSWUDy7fuZvnWPTz89iZaWgMXnBSkx3NeQdrxkBg1IJm4aN0H0RcpEESEgvQECtITuH5cYOLhg0daKNu2lxVVu1lRtYdlmxv4qzcOER1pjMxNPn4EcW5+KoUZiTrV1AfoslMR6ZDafYdYvnUPK6r2sKJqN6uq99J0pAWAhJhIRvTvx6gBKYwckMzI3GSG9++nI4kgpPsQRKTbtbQ6NtY1Ula9l9Xb97GmZh9rt++j8XAzAJERxpCsxEBI5CYfD4q0xBifKw9vug9BRLpdZIQxon8yI/on81mvrbXVUb37IGtqvJDYvo9Fm3bxwvJtx/cbkBLHyAHJFOf0Y1hOEsXZ/RianaSjiSCgQBCRbhMR8fcB62mjc4+379p/mDU1gYBYU7OPtTX7eHtDPUdbAmcozGBgegLF2YGQGJbTj+KcwKWzCoreo0AQkR6XkRTL5OIsJhdnHW872tLK5p0H2FC7nw21jWysa2Rj7X4Wrq+j2bvCKcJgUEYiQ7OTGOYFRGFmIkVZSaTER/vVnT5LgSAivoiOjKA4px/FOf34DH8/mjjS3MrmXQfYUNvIhtr9bKxtZENtI2+uqzt+KSwEbqoryko8HhCFmYkMyUqkID1B04N3kgJBRIJKTFQEw3L6MSyn3wntR5pbqdrdREX9ASp37qei/gAVOw/w5rp65pZWH98uwiA/LeF4WBRmJjIwPYFBGYnkpcZrNtjTUCCISEiIiYo4PiUH5Jywbt+ho1TWH6By5wEq6vdTsfMAFfUHWFLRwMGjLce3izAYkBrPoIxAQAxKT2BQRgID0xMZlJFAYmx4fySGd+9FpE9IjotmbEEqYwtST2h3zlHfeJgtDU1s2dXEll0HAq8NTcwvq2F309ETts9MimFQRuCIIj8tnoK0BPLTA6+5KXFERfbtowsFgoj0WWZGdnIc2clxjB+cfsr6vQePsnVXE1saAkFxbHlpZQMvrjhImyELIiOM3JQ4CtISKEiPJ997DbxPICspNuTv1lYgiEjYSomPZkx+CmPyU05Zd7SllZo9h6ja3URVQxPVuw8eX35rfT31jYdP2D4mMoIBqXHkpcWTlxpPXmrC8eX8tHj6p8QRHeRHGEETCGY2DbgfiATmOOfu9bkkEQlj0ZERx++paM+hoy3HQ6J690GqG5qo3nOQbbsPthsYEQY5yXHkHwuMtEBoDEgNtOWmxPs+hhEUgWBmkcCDwFSgGvjAzOY559b4W5mISPvioiMZmp3E0OykdtcfOtpCzd5DbNt9kG17mti2++DxwCjdspu/rao5fr/FMakJ0QxIiT9+ZDEgNY4BqceOOOLJ7OHTUkERCMAEoNw5VwFgZs8A0wEFgoiEpLjoyOOXvbanpdVRu+8Q2/ccZJv3s33PQbbvOcTWXU0s2rSL/d68UMfEREaQmxrHnZcP55qxA7q95mAJhDygqs37amDiyRuZ2SxgFsDAgQN7pzIRkR4QGWEMSI1nQGo87c40R+By2m27jwVF4Ahj+55DZPTQBIHBEgjtHQOdMg2rc242MBsCs532dFEiIn5KjosmOTeac3KTe+XvBcuQdzVQ0OZ9PrDdp1pERMJSsATCB0CxmRWaWQwwA5jnc00iImElKE4ZOeeazewbwKsELjt9zDm32ueyRETCSlAEAoBz7mXgZb/rEBEJV8FyykhERHymQBAREUCBICIiHgWCiIgAYM6F5v1dZlYPbAEygZ0+l+OncO5/OPcdwrv/4dx36Fr/BznnstpbEbKBcIyZlTrnPu7O7z4vnPsfzn2H8O5/OPcdeq7/OmUkIiKAAkFERDx9IRBm+12Az8K5/+Hcdwjv/odz36GH+h/yYwgiItI9+sIRgoiIdAMFgoiIACEUCGY2zczWm1m5md3dznozswe89avM7AI/6uwJHej7jV6fV5nZ+2Y21o86e8qZ+t9mu/Fm1mJm1/dmfT2pI303sylmtsLMVpvZ271dY0/qwP/9FDP7q5mt9Pr/FT/q7Alm9piZ1ZnZRx+zvvs/85xzQf9DYErsTUAREAOsBEaetM1VwHwCT1+bBCzxu+5e7PtFQJq3fGVf6XtH+99muzcJzJh7vd919+K/fSqBZ48P9N5n+113L/f/+8BPvOUsoAGI8bv2bur/pcAFwEcfs77bP/NC5QhhAlDunKtwzh0BngGmn7TNdOAJF7AYSDWz3N4utAecse/Oufedc7u9t4sJPHGur+jIvz3AN4E/A3W9WVwP60jfvwA875zbCuCcC7f+O6CfmRmQRCAQmukDnHPvEOjPx+n2z7xQCYQ8oKrN+2qv7Wy3CUVn26+bCXxr6CvO2H8zywOuA37Ti3X1ho782w8D0sxsoZktM7Obeq26nteR/v8aOIfAI3fLgDucc629U57vuv0zL2gekHMG1k7bydfLdmSbUNThfpnZJwkEwiU9WlHv6kj/fwnc5ZxrCXxR7DM60vcoYBxwGRAPLDKzxc65DT1dXC/oSP+vAFYAnwKGAK+b2bvOuX09XFsw6PbPvFAJhGqgoM37fALfCM52m1DUoX6Z2bnAHOBK59yuXqqtN3Sk/yXAM14YZAJXmVmzc+4vvVJhz+no//udzrkDwAEzewcYC/SFQOhI/78C3OsCJ9XLzawSGAEs7Z0SfdXtn3mhcsroA6DYzArNLAaYAcw7aZt5wE3eyPskYK9zrqa3C+0BZ+y7mQ0Enge+1Ee+GbZ1xv475wqdc4Odc4OB54Cv94EwgI79v38RmGxmUWaWAEwE1vZynT2lI/3fSuDoCDPLAYYDFb1apX+6/TMvJI4QnHPNZvYN4FUCVx485pxbbWZf89b/hsDVJVcB5UATgW8OIa+Dff8hkAE85H1LbnZ9ZCbIDva/T+pI351za83sFWAV0ArMcc61e5liqOngv/2PgN+bWRmBUyh3Oef6xLTYZvY0MAXINLNq4N+AaOi5zzxNXSEiIkDonDISEZEepkAQERFAgSAiIh4FgoiIAAoEERHxKBBERARQIIiIiOf/A7jfME3ZV37WAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate returns to obtain volatility\n",
    "returns = [(data['price'].iloc[i + 1] / data['price'].iloc[i]) - 1 for i in range(len(price) - 1)]\n",
    "\n",
    "# Check if volatility is consistent (assuming hourly volatility)\n",
    "sigma_24hour = calculate_std_dev(returns) * np.sqrt(24) #We multiply the std by str(24) because we got the volatility on one hour\n",
    "sigma_total=sigma_24hour*np.sqrt(252) #We multiply the std by str(252) because we got the volatility on one day\n",
    "print(f\"Annualized volatility: {sigma_total}\")\n",
    "\n",
    "# Calculate volume*signed to have the volume of the buying or selling order\n",
    "signed_volume = np.array([sign[i]*volume[i] for i in range(len(volume))])\n",
    "\n",
    "# Calculate price differences in steps of 2, we are using 2steps because the first price is before the transaction and the second one is after so we need this to have the impact of the order\n",
    "price_diff = np.array([price[i + 1] - price[i] for i in range(0, len(price) - 1, 2)])\n",
    "\n",
    "# Calculate squared of the volume*signed\n",
    "signed_volume_squared = np.array([sign[i] * volume[i] * volume[i] for i in range(len(price_diff))])\n",
    "\n",
    "# Define variables for the regression to find gamma and eta\n",
    "gamma_prices = price_diff\n",
    "eta_prices = -price_diff\n",
    "X_gamma = signed_volume.reshape(-1, 1)\n",
    "X_eta = np.column_stack((signed_volume, signed_volume_squared))\n",
    "\n",
    "# Create a linear regression model\n",
    "model = LinearRegression()\n",
    "\n",
    "# Regression for eta\n",
    "model.fit(X_eta, eta_prices)\n",
    "\n",
    "tau=1/24 #Tau is the number that represent at the step (we can do one transaction each hour)\n",
    "\n",
    "eta = model.coef_[1] * tau  # eta/tau is the coefficient for squared volumes, so multiply by tau to get eta\n",
    "\n",
    "# Predictions and evaluation for eta\n",
    "y_pred_eta = model.predict(X_eta)\n",
    "mse_eta = mean_squared_error(eta_prices, y_pred_eta)\n",
    "r2_eta = r2_score(eta_prices, y_pred_eta)\n",
    "\n",
    "print(f\"\\nEstimation of eta: {eta}\")\n",
    "print(f\"Mean Squared Error: {mse_eta} ; R^2: {r2_eta}\")\n",
    "\n",
    "# Regression for gamma\n",
    "model.fit(X_gamma, gamma_prices)\n",
    "gamma = model.coef_[0]  # Gamma is the slope of the regression line\n",
    "\n",
    "# Predictions and evaluation for gamma\n",
    "y_pred_gamma = model.predict(X_gamma)\n",
    "mse_gamma = mean_squared_error(gamma_prices, y_pred_gamma)\n",
    "r2_gamma = r2_score(gamma_prices, y_pred_gamma)\n",
    "\n",
    "print(f\"\\nEstimation of gamma: {gamma}\")\n",
    "print(f\"Mean Squared Error: {mse_gamma} ; R^2: {r2_gamma}\")\n",
    "\n",
    "# Liquidation strategy\n",
    "X_to_liquidate=15000 #It's an example\n",
    "possible_transaction=[(1/24)*i for i in range(1,25)] #One transaction possible each hour\n",
    "Lambda_value = 10**(-4.0) #It's an example of the aversion (taste of risk)\n",
    "K_value = np.sqrt(Lambda_value * (sigma_total**2) / eta)\n",
    "liquidation_quantity=[]\n",
    "for t_value in possible_transaction :\n",
    "    liquidation_quantity.append(np.sinh(K_value * (1 - t_value)) * X_to_liquidate / np.sinh(K_value * 1))\n",
    "\n",
    "# Plot the liquidation strategy\n",
    "plt.plot(possible_transaction, liquidation_quantity)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question E (Q4 of TD5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Importation and transformation fo the data \n",
    "data=pd.read_excel(\"Dataset TD5.xlsx\")\n",
    "data_GBP=data.iloc[2:, :3]\n",
    "data_SEK=data.iloc[2:, 4:7]\n",
    "data_CAD=data.iloc[2:, 8:11]\n",
    "\n",
    "#We need to have the half of the bid-ask spread on each FX rates to get the returns of each one.\n",
    "half_bidask_GBP=[(data_GBP.iloc[i,1]+data_GBP.iloc[i,2])/2 for i in range(0,len(data_GBP)-1)]\n",
    "half_bidask_SEK=[(data_SEK.iloc[i,1]+data_SEK.iloc[i,2])/2 for i in range(0,len(data_SEK)-1)]\n",
    "half_bidask_CAD=[(data_CAD.iloc[i,1]+data_CAD.iloc[i,2])/2 for i in range(0,len(data_CAD)-1)]\n",
    "\n",
    "#Returns of each FX rates\n",
    "returns_gbp=[(half_bidask_GBP[i+1]/half_bidask_GBP[i])-1 for i in range(0,len(data_GBP)-2)]\n",
    "returns_sek=[(half_bidask_SEK[i+1]/half_bidask_SEK[i])-1 for i in range(0,len(data_GBP)-2)]\n",
    "returns_cad=[(half_bidask_CAD[i+1]/half_bidask_CAD[i])-1 for i in range(0,len(data_GBP)-2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We tried to do the correlation matrix by using the wavelets but we didn't get a good result, this is the part on the wavelets:\n",
    "def calculate_empirical_moment(data, k, N):\n",
    "    result = 0\n",
    "    for i in range(1, len(data)*N):\n",
    "        result += np.abs(data[int(i/N)] - data[int((i-1)/N)])**k\n",
    "    return result * (1/N * len(data))\n",
    "\n",
    "def calculate_empirical_moment_halved_resolution(data, k, N):\n",
    "    result = 0\n",
    "    for i in range(1, int((len(data)*N)/2)):\n",
    "        result += np.abs(data[int(2*i/N)] - data[int(2*(i-1)/N)])**k\n",
    "    return result * (2/N * len(data))\n",
    "\n",
    "#We didn't get any good result from this so we did the correlation matrix with excel and the result is : \n",
    "correlation_matrix = np.array([[1, 0.191420071, 0.244926851],\n",
    "                               [0.191420071, 1, 0.192442093],\n",
    "                               [0.244926851, 0.192442093, 1]]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5\n",
      "0.5\n",
      "0.5\n",
      "0.0008814409953038501\n",
      "0.00046263269044024267\n",
      "0.0007159859130200779\n",
      "Covariance Matrix:\n",
      "[[7.76938228e-07 7.80579310e-08 1.54573173e-07]\n",
      " [7.80579310e-08 2.14029006e-07 6.37442282e-08]\n",
      " [1.54573173e-07 6.37442282e-08 5.12635828e-07]]\n"
     ]
    }
   ],
   "source": [
    "#This is the part on the Hursr exhibitor \n",
    "\n",
    "#Method to have the moment of the data for each k and N each one > 0 \n",
    "def calculate_empirical_moment(data, k, N):\n",
    "    result = 0\n",
    "    for i in range(1, len(data)*N):\n",
    "        result += np.abs(data[int(i/N)] - data[int((i-1)/N)])**k\n",
    "    return result * (1/N * len(data))\n",
    "\n",
    "#Method to have the moment halved resolution for k = 2 and any N > 0 \n",
    "def calculate_empirical_moment_halved_resolution(data, k, N):\n",
    "    result = 0\n",
    "    for i in range(1, int((len(data)*N)/2)):\n",
    "        result += np.abs(data[int(2*i/N)] - data[int(2*(i-1)/N)])**k\n",
    "    return result * (2/N * len(data))\n",
    "\n",
    "#Method to estimate the Hurst exhibitor\n",
    "def estimate_Hurst(data, k, N):\n",
    "    M2 = calculate_empirical_moment(data, k, N)\n",
    "    M2_halved = calculate_empirical_moment_halved_resolution(data, k, N)\n",
    "    return 0.5 * math.log2(M2_halved/M2)\n",
    "\n",
    "#Method to get the scale volatily by using the Hurst exhibitor\n",
    "def scale_volatility(std, hurst_exponent):\n",
    "    return std * (2 ** hurst_exponent)\n",
    "\n",
    "hurst_gbp= estimate_Hurst(returns_gbp, 2, 10)\n",
    "hurst_sek=estimate_Hurst(returns_sek, 2, 10)\n",
    "hurst_cad=estimate_Hurst(returns_cad, 2, 10)\n",
    "print(hurst_gbp)\n",
    "print(hurst_sek)\n",
    "print(hurst_cad)\n",
    "scale_vol=[]\n",
    "scale_vol.append(scale_volatility(calculate_std_dev(returns_gbp),hurst_gbp))\n",
    "scale_vol.append(scale_volatility(calculate_std_dev(returns_sek),hurst_sek))\n",
    "scale_vol.append(scale_volatility(calculate_std_dev(returns_cad),hurst_cad))\n",
    "for i in scale_vol:\n",
    "    print(i)\n",
    "\n",
    "#Method to get the mutiplication between the correlation_matrix obtained with Excel and the matrix of the scale volatility\n",
    "def calculate_covariance_matrix(correlation_matrix, volatility_vector):\n",
    "    correlation_matrix = np.array(correlation_matrix)\n",
    "    Sigma = np.diag(volatility_vector)\n",
    "    covariance_matrix = np.dot(Sigma, np.dot(correlation_matrix, Sigma))\n",
    "    return covariance_matrix\n",
    "covariance_matrix = calculate_covariance_matrix(correlation_matrix, scale_vol)\n",
    "print(\"Covariance Matrix:\")\n",
    "print(covariance_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#To calculate the volatility directly from the series of price we decided to try to create different datasets for different time scales\n",
    "#Create windows (5)\n",
    "column_names = ['Date', 'High', 'Low']\n",
    "\n",
    "data_GBP.columns = column_names\n",
    "data_SEK.columns = column_names\n",
    "data_CAD.columns = column_names\n",
    "\n",
    "data_GBP['Date'] = pd.to_datetime(data_GBP['Date'])\n",
    "data_SEK['Date'] = pd.to_datetime(data_SEK['Date'])\n",
    "data_CAD['Date'] = pd.to_datetime(data_CAD['Date'])\n",
    "\n",
    "data_GBP.index = pd.to_datetime(data_GBP.index)\n",
    "data_SEK.index = pd.to_datetime(data_SEK.index)\n",
    "data_CAD.index = pd.to_datetime(data_CAD.index)\n",
    "\n",
    "date_ranges = [\n",
    "    ('2016-03-07 09:00:00', '2016-03-07 10:00:00'),\n",
    "    ('2016-03-07 09:00:00', '2016-03-07 19:00:00'),\n",
    "    ('2016-03-07 09:00:00', '2016-03-08 09:00:00'),\n",
    "    ('2016-03-07 09:00:00', '2016-03-10 09:00:00'),\n",
    "    ('2016-03-07 09:00:00', '2016-03-14 09:00:00'),\n",
    "    ('2016-03-07 09:00:00', '2016-04-07 09:00:00'),\n",
    "    ('2016-03-07 09:00:00', '2016-06-07 09:00:00')\n",
    "]\n",
    "#After create this samples, we didn't succeed to create the dataset with the final date that are each values on the right\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result with Overlapping : \n",
      "\n",
      "- Standard deviation of GBP : 0.0006232729049951721\n",
      "- Standard deviation of GBP annualised: 0.002159080676865376\n",
      "- Standard deviation of SEK: 0.00032713071260887243\n",
      "- Standard deviation of SEK annualised: 0.0011332140299095596\n",
      "- Standard deviation of CAD: 0.00032713071260887243\n",
      "- Standard deviation of CAD annualised: 0.0011332140299095596\n",
      "\n",
      "\n",
      "Result with Non-Overlapping : \n",
      "\n",
      "- Standard deviation of GBP : 0.03799434923003744\n",
      "- Standard deviation of GBP annualised: 0.13161628653388058\n",
      "- Standard deviation of SEK : 0.012694382028994778\n",
      "- Standard deviation of SEK annualised: 0.04397462928981649\n",
      "- Standard deviation of CAD : 0.010946213754529869\n",
      "- Standard deviation of CAD annualised: 0.037918796746710025\n"
     ]
    }
   ],
   "source": [
    "#We decided to just apply the two method (overlapping and non-overlapping) on the whole dataset\n",
    "\n",
    "#With overlapping\n",
    "std_GBP = calculate_std_dev(returns_gbp)\n",
    "std_SEK = calculate_std_dev(returns_sek)\n",
    "std_CAD = calculate_std_dev(returns_sek)\n",
    "print(\"Result with Overlapping : \\n\")\n",
    "print(f\"- Standard deviation of GBP : {std_GBP}\")\n",
    "print(f\"- Standard deviation of GBP annualised: {std_GBP*np.sqrt(12)}\") #Annualized (we have 6 months so for annualised it's the root of 12)\n",
    "print(f\"- Standard deviation of SEK: {std_SEK}\")\n",
    "print(f\"- Standard deviation of SEK annualised: {std_SEK*np.sqrt(12)}\") #Annualised\n",
    "print(f\"- Standard deviation of CAD: {std_CAD}\")\n",
    "print(f\"- Standard deviation of CAD annualised: {std_CAD*np.sqrt(12)}\") #Annualised\n",
    "print()\n",
    "\n",
    "#Without overlapping\n",
    "returns_gbp_without_overlapping=[(half_bidask_GBP[i+1]/half_bidask_GBP[0])-1 for i in range(0,len(data_GBP)-2)]\n",
    "returns_sek_without_overlapping=[(half_bidask_SEK[i+1]/half_bidask_SEK[0])-1 for i in range(0,len(data_GBP)-2)]\n",
    "returns_cad_without_overlapping=[(half_bidask_CAD[i+1]/half_bidask_CAD[0])-1 for i in range(0,len(data_GBP)-2)]\n",
    "std_GBP_without_overlapping = calculate_std_dev(returns_gbp_without_overlapping)\n",
    "std_SEK_without_overlapping = calculate_std_dev(returns_sek_without_overlapping)\n",
    "std_CAD_without_overlapping = calculate_std_dev(returns_cad_without_overlapping)\n",
    "print(\"Result with Non-Overlapping : \\n\")\n",
    "print(f\"- Standard deviation of GBP : {std_GBP_without_overlapping}\")\n",
    "print(f\"- Standard deviation of GBP annualised: {std_GBP_without_overlapping*np.sqrt(12)}\") #Annualised\n",
    "print(f\"- Standard deviation of SEK : {std_SEK_without_overlapping}\")\n",
    "print(f\"- Standard deviation of SEK annualised: {std_SEK_without_overlapping*np.sqrt(12)}\") #Annualised\n",
    "print(f\"- Standard deviation of CAD : {std_CAD_without_overlapping}\")\n",
    "print(f\"- Standard deviation of CAD annualised: {std_CAD_without_overlapping*np.sqrt(12)}\") #Annualised"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}