14_Value_at_Risk.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "7ddf853b-edf6-4db4-8e72-abbf1657cccd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import yfinance as yf\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "from tqdm.notebook import tqdm\n",
    "from scipy.stats import shapiro, normaltest, kstest\n",
    "\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "541bbe79-2d14-4bcd-a3c6-138b6ad7e9bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "stocks = [\n",
    "    \"PGSUS.IS\",\n",
    "    \"ISCTR.IS\",\n",
    "    \"KCHOL.IS\",\n",
    "    \"ASELS.IS\",\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "09d000f2-601d-48ac-9d87-4e7270707bfb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[*********************100%***********************]  4 of 4 completed\n"
     ]
    }
   ],
   "source": [
    "data = yf.download(tickers=stocks, period=\"max\",\n",
    "                   start=\"2020-01-03\")['Adj Close']\n",
    "\n",
    "data.clip(lower=0, inplace=True)\n",
    "data.fillna(0, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "b86c37a8-bd4e-42e8-ae5a-94bb49b1bb5b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ASELS.IS</th>\n",
       "      <th>ISCTR.IS</th>\n",
       "      <th>KCHOL.IS</th>\n",
       "      <th>PGSUS.IS</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2021-08-31</th>\n",
       "      <td>0.005060</td>\n",
       "      <td>0.019504</td>\n",
       "      <td>-0.023729</td>\n",
       "      <td>-0.000708</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2021-09-01</th>\n",
       "      <td>-0.011328</td>\n",
       "      <td>-0.005217</td>\n",
       "      <td>0.007812</td>\n",
       "      <td>0.012748</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2021-09-02</th>\n",
       "      <td>0.001273</td>\n",
       "      <td>0.005245</td>\n",
       "      <td>-0.001723</td>\n",
       "      <td>-0.002797</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2021-09-03</th>\n",
       "      <td>-0.001907</td>\n",
       "      <td>-0.010435</td>\n",
       "      <td>-0.006040</td>\n",
       "      <td>-0.005610</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2021-09-06</th>\n",
       "      <td>-0.000637</td>\n",
       "      <td>0.008787</td>\n",
       "      <td>0.004340</td>\n",
       "      <td>0.005642</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            ASELS.IS  ISCTR.IS  KCHOL.IS  PGSUS.IS\n",
       "Date                                              \n",
       "2021-08-31  0.005060  0.019504 -0.023729 -0.000708\n",
       "2021-09-01 -0.011328 -0.005217  0.007812  0.012748\n",
       "2021-09-02  0.001273  0.005245 -0.001723 -0.002797\n",
       "2021-09-03 -0.001907 -0.010435 -0.006040 -0.005610\n",
       "2021-09-06 -0.000637  0.008787  0.004340  0.005642"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "returns = data.pct_change().dropna()\n",
    "returns.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "f17639b0-461f-4e16-aa3e-cdd95f27fd35",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Portfolio with random weights\n",
    "weights = np.random.random(len(stocks))\n",
    "weights = weights / np.sum(weights)\n",
    "\n",
    "# Portfolio PnLs\n",
    "PnL = (weights * returns).sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "b5d5f7dd-e20f-4d04-9e0a-b8fc5e9438d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Confidence Level    Value-at-Risk\n",
      "---------------------------------\n",
      "            99.9          0.10379\n",
      "            99.0          0.05813\n",
      "            95.0          0.03256\n"
     ]
    }
   ],
   "source": [
    "percs = [0.1, 1, 5]\n",
    "var = np.percentile(PnL, percs, interpolation=\"lower\")\n",
    "print('%16s %16s' % ('Confidence Level', 'Value-at-Risk'))\n",
    "print(33 * '-')\n",
    "for pair in zip(percs, var):\n",
    "    print('%16.1f %16.5f' % (100 - pair[0], -pair[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "0b7cca1c-f135-4e89-9af6-890e88f2bb4f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Historical 99% VaR: -0.05813094911174911\n",
      "Historical 95% VaR: -0.03255777292807174\n"
     ]
    }
   ],
   "source": [
    "# Historical VaRs\n",
    "# https://medium.com/@matt_84072/historical-simulation-value-at-risk-explained-with-python-code-a904d848d146\n",
    "var_1p = np.percentile(PnL, 1, interpolation=\"lower\")\n",
    "var_5p = np.percentile(PnL, 5, interpolation=\"lower\")\n",
    "\n",
    "print(f\"Historical 99% VaR: {var_1p}\")\n",
    "print(f\"Historical 95% VaR: {var_5p}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "f9c39e4a-d26f-4c8e-bed8-0998949d1533",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x238214c12e0>"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x2382188a340>"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x238214afa30>"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "[Text(0.5, 6.79999999999999, 'PnL'),\n",
       " Text(0.42499999999999716, 0.5, 'Probability')]"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAFuCAYAAABQqakCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABDFElEQVR4nO3deXxU9b3/8dcn+76SBEjCvgi4oCKCdcFavWBdqqJCtUqrP7TW3nJbf1Z7q7daf/XaW1vbW3utt25VW7RWrAsudaFWQQURF4RAWJNAJivZ9/n+/pgTHEKWSZgzZ87M5/l45MHMWd+ZhA+H7/d7vkeMMSillHJGjNMBlFIqmmkRVkopB2kRVkopB2kRVkopB2kRVkopB8U5HSBYFi5caF555RWnY7jbggW+P9escTJFvxY8ugCANcvWOJpDqRGSgVZEzJVwTU2N0xGUUmrYIqYIK6WUG2kRVkopB0VMm7CKbNoWrCKVFmGlVEC6urooLy+nvb3d6ShhKykpiaKiIuLj4wPeR4uwcoVfrP0FADedcpPDSaJXeXk56enpTJgwAZEBO/ujljGG2tpaysvLmThxYsD7aZuwcoUXt73Ii9tedDpGVGtvbyc3N1cL8ABEhNzc3GH/T0GLsFIqYFqABzeSz0eLsFJKOUiLsFLKNb71rW+Rn5/P0UcffcjyH/7whxx77LFcddVVB5c98cQT3Hffff0eZ9KkSZSUlByybMWKFdxzzz0DnnvChAkcc8wxHHvssZxxxhns2bNn5N+IHy3CyhWS45NJjk92OoZy2LJly+g7PUFDQwMbN27kk08+ISEhgU8//ZS2tjYeeeQRvvOd7/R7nCVLlrBy5cqD771eL8888wxLliwZ9PxvvfUWn3zyCQsWLOCuu+468m8ILcLKJV6+4mVevuJlp2Moh51++unk5OQcsiwmJoauri6MMbS2thIfH88vfvELvvvd7w44VGzp0qU89dRTB9+//fbbjB8/nvHjx/O1r32NE088kVmzZvHggw/2u//8+fOpqKgIyvekQ9SUUsO3YgVs2hTcY86eDQM0HwwmPT2dc889l+OPP56zzjqLzMxM3n//fW677bYB9znmmGOIiYnh448/5rjjjmPlypUsXboUgIcffpicnBza2to46aSTuOSSS8jNzT1k/1deeYWvfe1rw87aH70SVq7w03/8lJ/+46dOx1Bh6uabb2bTpk3ce++93Hbbbdx555384Q9/4LLLLhuw2WDp0qWsXLmS7u5unnvuOS699FIAfvOb33Dccccxb948ysrK2L59+8F9zjzzTAoLC3n55ZcPFu0jpVfCyhXe2PUGALedMfDVjQqhEVyxhsJHH32EMYbp06dz66238uqrr/LNb36T7du3M3Xq1EO2XbJkCeeccw5nnHEGxx57LAUFBaxZs4bXX3+ddevWkZKSwoIFCw4Z9/vWW2+RlZXFFVdcwX/8x3/wy1/+8ogzaxFWKgROPuVUKj2eAdePLijg/bXvhDBRZLrtttt48MEH6erqoqenB/C1Gbe2th627eTJkxk1ahS33HIL3/ve9wBfJ192djYpKSls3bqV995777D94uLiuO+++zjmmGP48Y9/fFgb9XBpEVYqBCo9HlY88NKA6++7/qshTONeS5cuZc2aNdTU1FBUVMQdd9zBNddcA8Bzzz3HnDlzGDt2LACzZ88+OKTsuOOOG/B4t9xyCxdffDEACxcu5IEHHmDGjBlMnz6defPm9bvfmDFjWLp0Kffff/+gbc+BEGPMER0gXMyZM8ds2LDB6Rjupk/WsM34yVOHLMJ7dmwfcH042LJlCzNmzHA6Rtgb4HMa8FY6vRJWrpCbkjv0Rkq5kBZh5Qp/veyvTkdQyhY6RE0ppRykRVi5wq2v38qtr9/qdAylgs7WIiwiC0WkRERKReSWftafLiIbRaRbRBb3WTdORF4TkS0i8rmITLAzqwpv68rXsa58ndMxlAo624qwiMQC9wOLgJnAUhGZ2WezvcAy4E/9HOKPwH8ZY2YAc4Equ7IqpZRT7LwSnguUGmN2GmM6gZXAhf4bGGN2G2M+Abz+y61iHWeM+bu1XbMx5vDR1kqpqNI7neTs2bOZM2fOweU6lWX/CoEyv/fl1rJATAMOiMizIvKRiPyXdWV9CBFZLiIbRGRDdXV1ECIrpcLdW2+9xaZNm+i9L0CnsrRHHHAacBNwEjAJX7PFIYwxDxpj5hhj5uTl5YU2oQqpoowiijKKnI6hwpBOZTmwCqDY732RtSwQ5cAmY8xOABF5DpgHPBTMgMo9nrj4CacjqD5672L0d9msy7jhpBto7Wrl3CfPPWz9stnLWDZ7GTWtNSx++pC++IDuhhQRzjnnHESE6667juXLl7t+Kks7i/B6YKqITMRXfJcAXx/GvlkikmeMqQa+DOg9yUpFuXfeeYfCwkKqqqo4++yzOeqoozj99NO5+eabufnmmwG49tprD05l+dprr3Hsscfy4x//+LBj9U5lOWvWLJ577jnuuOMOwDeV5apVqwAOTmXZW4TPPPNM6urqSEtL46c/Dc7UqrYVYWNMt4jcCLwKxAIPG2M2i8idwAZjzPMichKwCsgGzheRO4wxs4wxPSJyE/CG+B5f+iHwv3ZlVeFvxSsrALhv4X2O5lBfGOzKNSU+ZdD1o1JGjWgekMJCX7dSfn4+F110ER988AGnn376wfU6lWUfxpjVwOo+y273e70eXzNFf/v+HTjWznzKPTZVbnI6gnJYS0sLXq+X9PR0WlpaeO2117j99tsP2caNU1mGa8ecUkodwuPxcOqpp3Lccccxd+5cvvrVr7Jw4cKD6/2nsszKyjo4lWV7e/ugU1lu3br1kKksu7u7mTFjBrfccktAU1keKZ3KUn1Bp7IcsaEmbfd4PNy9auDfT53KMnLoVJZKOWCoSdtvvuCEEKZRbqJFWLnCtNxpTkdQyhZahJUrPHh+/4PmVWgZY/ANWFL9GUnzrnbMKaUCkpSURG1t7YgKTTQwxlBbW0tSUtKw9tMrYeUKy19YDugVsZOKioooLy9H52kZWFJSEkVFw7u9XouwcoVttducjhD14uPjmThxotMxIo42RyillIO0CCullIO0CCullIO0TVi5wuzRs52OoJQttAgrV9DZ01Sk0uYIpZRykBZh5QpXPnslVz57pdMxlAo6bY5QrlDeWO50BKVsoVfCSinlIC3CSinlIC3CSinlIG0TVq4wv2i+0xGUsoUWYeUKd3/lbqcjKGULbY5QSikHaRFWrnDJ05dwydOXOB1DqaDT5gjlCrWttU5HUMoWeiWslFIO0iKslFIO0iKslFIO0jZh5QpnTTzL6QhK2cLWK2ERWSgiJSJSKiK39LP+dBHZKCLdIrK4n/UZIlIuIr+1M6cKf7edcRu3nXGb0zGUCjrbirCIxAL3A4uAmcBSEZnZZ7O9wDLgTwMc5qfA23ZlVEopp9l5JTwXKDXG7DTGdAIrgQv9NzDG7DbGfAJ4++4sIicCBcBrNmZULrHoyUUsenKR0zGUCjo7i3AhUOb3vtxaNiQRiQHuBW4aYrvlIrJBRDZUV1ePOKgKf21dbbR1tTkdQ6mgC9fRETcAq40xg87kbYx50BgzxxgzJy8vL0TRlFIqeOwcHVEBFPu9L7KWBWI+cJqI3ACkAQki0myMOaxzTyml3MzOIrwemCoiE/EV3yXA1wPZ0RhzRe9rEVkGzNECrJSKRLYVYWNMt4jcCLwKxAIPG2M2i8idwAZjzPMichKwCsgGzheRO4wxs+zKpNzrvGnnOR1BKVvYerOGMWY1sLrPstv9Xq/H10wx2DEeBR61IZ5ykZtOGbSPNmy1d/WwYU896XMvob6lk+zUBKcjqTATrh1zSrlea2c3f91YzsY99aTPW8xfPiynqb3L6VgqzGgRVq6w4NEFLHh0gdMxhuXtbTXUt3Zx4eyxVD1+Ez1ew0uf7scY43Q0FUa0CCtlg/rWTrZ5mphdlMX43FS66ys4Y3oensYOdtW2OB1PhREtwkrZYMPuemJihOPHZR1cNr0gnbTEODaVHXAslwo/WoSVCrKuHi8lniZmjE4nNfGLvu/YGOHYokzK6tqobe5wMKEKJ1qElQqy3bUt9HgN0wrSD1s3a2wGIlDiaXIgmQpHOp+wcoXLZl3mdISA7ahqISk+hsKs5MPWpSTEMTYzmZ01LZwyeZQD6VS40SKsXOGGk25wOkJAur1edtW0MCU/jZgY6XebSXmp/HN7DY1tXWQkx4c4oQo32hyhXKG1q5XWrlanYwxp34F2Onu8TM5PHXCbiaN863bW6CgJpUVYucS5T57LuU+e63SMIVXUtyHQb1NEr+yUBLJT4tmlRVihRVipoKo40EZeeiKJcbGDbleck8L+hja8Xr1xI9ppEVYqSLp7vFQ2tlOYPfBVcK/CrGS6egxVOlQt6mkRVipIPI0d9HgNRYM0RfTqba7Yd0CfFhLttAgrFSQVVkEdG0ARTk2MIzM5nop6LcLRToeoKVdYNnuZ0xGGVNnYTk5KAknxg7cH9yrMSmZndbNO6BPltAgrV3BDEa5qbKc4JyXg7cdmJfH5/kYOtOr0ltFMmyOUK9S01lDTWuN0jIElZdDS2UN+emLAuxRkJAHgaWq3K5VyAb0SVq6w+OnFAKxZtsbZIAPJ9j3TtrewBiInJYG4GMHTqCMkopleCSsVDNnFCJA3jCvhmBghLz2Rqka9Eo5mWoSVCobsInJSE4iPHd5fqYL0JKqaOoD+55lQkU+LsFLBkF08rPbgXgUZiXR7DWQU2BBKuYEWYaWOUF1LJ5KcyagRFOH83jbk7HFBTqXcQjvmlCt8e863nY4woJJK3wTtuSN4nH1WSjxxMUJX1thgx1IuoUVYucLlR1/udIQBbbOekpGbNvwr4RgRclIT8GSMDnYs5RLaHKFcoayhjLKGMqdj9GtrZROmo4XUhMDulOsrNy0BMscEOZVyCy3CyhW+seobfGPVN5yO0a9tniZo3I/IyEY4jEpNRJIzqWvpDHIy5QZahJU6AsYYtlU2QcP+ER8jN83XltzbtqyiixZhpY7AvoZ2mjq6j6gIj7LakksqG4MVS7mIrUVYRBaKSImIlIrILf2sP11ENopIt4gs9ls+W0TWichmEflERMK3V0ZFtd5OuSMpwikJsZiOFko8eiUcjWwrwiISC9wPLAJmAktFZGafzfYCy4A/9VneClxljJkFLATuE5Esu7IqNVI7q63nxDVVjfgYIgIN+7U5IkrZOURtLlBqjNkJICIrgQuBz3s3MMbsttZ5/Xc0xmzze71PRKqAPOCAjXlVGPvB/B84HaFfO6ubyUyO50DnET60s3E/2zy+uYVH2sGn3MnOIlwI+I8pKgdOHu5BRGQukADs6GfdcmA5wLhxesdRJDt/+vlOR+jXrpoWJo5K5aMjPVDDfpo7uqk40EZRduBzEiv3C+uOOREZAzwOfNMY4+273hjzoDFmjjFmTl5eXugDqpApqSmhpKbE6RiH2VndwqS81CM/UMM+QEdIRCM7i3AFUOz3vshaFhARyQBeAv7dGPNekLMpl7nuxeu47sXrnI5xiJaObiob25k0KhhFuBLw3fihooudRXg9MFVEJopIArAEeD6QHa3tVwF/NMY8Y2NGpUZsd62vHXhSXtqRH6y7ncKs5C9GW6ioYVsRNsZ0AzcCrwJbgKeNMZtF5E4RuQBARE4SkXLgUuD3IrLZ2v0y4HRgmYhssr5m25VVqZHoHRkxMRhXwsC0gjRtjohCtk7gY4xZDazus+x2v9fr8TVT9N3vCeAJO7MpdaR21QS3CE8fncE7pTV09XiHPTm8ci/9SSs1QjurmynMSg74EfdDmT46ja4ew+6aIxzuplxFp7JUrvDj03/sdITD7KoJ0sgIy5S8dABKq5qZWpAetOOq8KZFWLnCVyZ9xekIhzDGsLO6hYtOKAzaMSfn+wr69qpmFgXtqCrcaXOEcoVNlZvYVLnJ6RgH1TR30tTRHZzhaZaUhDgKs5IprWoO2jFV+NMrYeUKK15ZAcCaZWsczdFrZ7WvUE4MxvA0P1Py07QIRxm9ElZqBHpHRgTzShhgan4aO6qb6fGaoB5XhS8twkqNwM6aFhLiYhiblRzU407JT6Oj20tFfVtQj6vClxZhpUZgZ3ULE3NTiY0J7oxnU/J9zRul1XrTRrTQIqzUCOysaQ7aTRr+eovwdo+2C0cL7ZhTrvCzs37mdISDeryGsrpWzpkZ/MfUZ6UkMCotUTvnoogWYeUKpxSf4nSEg/Y3tNHVYxifa8+8v1PyUymt1iIcLbQ5QrnC2rK1rC1b63QMAMrqfJ1m43LsKcJT89MprfI9ZUNFPr0SVq7wozd+BITHOOGyulbAviI8JT+NpvZuqpo6KMhIsuUcKnzolbBSw7S3rpXYGGFMpj0F8uAICW0XjgpahJUapr11rRRmJRNn03STWoSjixZhpYZpb12rbU0RAPnpiaQnxbG9SscKRwMtwkoNU1ldK8U5wb1Tzp+I6BwSUUQ75pQr3LfwPqcjAL6He9a2dFJs45UwwJS8NN4qqbb1HCo86JWwcoXZo2cze/Rsp2NQVm/vyIheUwvSqGnu4EBrp63nUc7TIqxc4fWdr/P6ztedjsHe2tAUYe2cix7aHKFc4a637wKcf8LGXpvHCPfyf9TRnAk5tp5LOUuvhJUahrK6VtKT4shMjrf1PIXZySTFx+iVcBTQIqzUMPQOTxMJ7hSWfcXGCJNGpbFdi3DE0yKs1DDsrWulONvepoheOkwtOmgRVipAXq+hvL6NcTbNntbX1Pw0Kg600drZHZLzKWdox5xyhd+f93unI1Dd3EFHt9f2McK9ekdI7Khq4ZiizJCcU4WeFmHlCtNHTXc6QshGRvTyf9SRFuHIpc0RyhVeKHmBF0pecDRDqMYI9xqfm0pcjGi7cISztQiLyEIRKRGRUhG5pZ/1p4vIRhHpFpHFfdZdLSLbra+r7cypwt+96+7l3nX3Opphb10rIlAY5CcsDyQhLobxuSn6vLkIZ1tzhIjEAvcDZwPlwHoRed4Y87nfZnuBZcBNffbNAf4DmAMY4ENr33q78io1mJNPOZXKcV+B/ClMnX5404jH4zmi43sqPYyfPPXwFfO/SWnGaMbfsZTRBQW8v/adIzqPCj92tgnPBUqNMTsBRGQlcCFwsAgbY3Zb67x99v0X4O/GmDpr/d+BhcCfbcyr1IAqPR4Kzz8NEWHxAy8dtv7mC044ouP3eL2s6Oe4a3fUsGFPPd/53Yv89w3nHdE5VHgKqDlCRM4XkeE2XRQCZX7vy61lQdtXRJaLyAYR2VBdrTNOKXs1tHfZfqdcXzmpCRiDTuQTwQItrJcD20Xk5yJylJ2BhsMY86AxZo4xZk5eXp7TcVQki4mnpaMn9EU4JQGAuhYtwpEqoOYIY8yVIpIBLAUeFREDPAL82Rgz0PT/FUCx3/sia1kgKoAFffZdE+C+KgI9ftHjzgZI9U2ik5Ec2lGd2alWEdYr4YgVcBODMaYReAZYCYwBLgI2ish3B9hlPTBVRCaKSAKwBHg+wNO9CpwjItkikg2cYy1TUao4s5jizOKhN7RLai5AyK+E42NjyEiK0yvhCBZom/CFIrIK39VoPDDXGLMIOA74QX/7GGO6gRvxFc8twNPGmM0icqeIXGAd9yQRKQcuBX4vIputfeuAn+Ir5OuBO3s76VR0euqzp3jqs6ecC5DmTBEG39WwFuHIFej/rS4GfmWMedt/oTGmVUSuGWgnY8xqYHWfZbf7vV6Pr6mhv30fBh4OMJ+KcP+z4X8AuPzoy50JkJpLXIyQHB8b8lPnpCZQXt8G2Dtzm3JGoM0RlX0LsIjcA2CMeSPoqZQKN6m5ZCbH2z6FZX9yUhPo8ZqD7dIqsgRahM/uZ9miYAZRKqxZRdgJvSMkyChw5PzKXoM2R4jIt4EbgMki8onfqnTgXTuDKRUujDGQlkuGU0XYGiFBuhbhSDRUm/CfgJeBuwH/uR+atKNMRYvalk4kLtGxK+Gk+FhSEmJpyRjtyPmVvYYqwsYYs1tEvtN3hYjkaCFWofLMZc84du7eKSydKsLguxpu0eaIiBTIlfB5wIf4JtLx75UwwCSbcil1iFEpoxw7d5lVhDOSnJt+OyclgbKMAowxjnQOKvsM+ltljDnP+nNiaOIo1b9HNz0KwLLZy0J+7oNF2OErYYlPpqqpg4KMJMdyqOAbqmNu0KmhjDEbgxtHqf45WYT31rVi2hqIj3XuGQi9nXOlVc1ahCPMUP+/GmwWbQN8OYhZlApLe+taoaXW0Qy9RXi7p4kvTXGuaUYF31DNEWeGKohS4aqsrg2anS3CKQmxmM42Sqv1KRuRZqjmiC8bY94UkYv7W2+MedaeWEqFh85uL/sa2hy/EhYRaKyktCrQKbmVWwzVHHEG8CZwfj/rDKBFWEW0igNtGIPjRRiAJo8+9DMCDdUc8R/Wn98MTRyl+rf6itVDb2SD3pERYVGEGz3UNHdyoLWTrN5bmZXrBTqVZa6I/MZ6MvKHIvJrEcm1O5xSvVLiU0iJD82j5v313qjhdJswAI2VAHo1HGECHXOzEqgGLgEWW68dnNxVRZvfrf8dv1v/u5Cft6yulYS4GGhvDPm5D9Poe6Lzdi3CESXQIjzGGPNTY8wu6+suQO+hVCHz9OaneXrz0yE/7966Voqzk/F1gTistZ6k+Bi9Eo4wgRbh10RkiYjEWF+XoY8bUlFgb10r43JC3wzSP8PkvDQtwhFmqCFqTXwxZ8QK4AlrVQzQDNxkZzilnLa3rpU547OdjnHQlPw0NuyudzqGCqJBr4SNMenGmAzrzxhjTJz1FWOMyQhVSKWc0NDaRVN7N8VhcyUMU/LSqDjQRktHt9NRVJAEPC2U9dTjqcDBG9f7PvJIqUjSOzIinIrw1II0AHZWt3BMUabDaVQwBFSEReRa4Hv4Hsq5CZgHrEPnjlAhsmbZmpCfs7cIh0+bsK85AqC0ukmLcIQItGPue8BJwB5rPonjgQN2hVIqHITjlfD43FTiYoRtHu2cixSBFuF2Y0w7gIgkGmO2AtPti6XUoX6x9hf8Yu0vQnrOvXWt5KYmkJbo3GTufcXHxjAlP40t+8Ng3LIKikCLcLmIZAHPAX8Xkb8Be+wKpVRfL257kRe3vRjSc5bVtVIURlfBvWaOzWDzPi3CkSKgImyMucgYc8AY8xPgNuAh4Gs25lLKcXvqWpiQG35FeNbYTKqbOqhqanc6igqCgB8VICIniMi/AscC5caYTvtiKeWszm4vFfVtjA/DK+FZY32jQz/Xq+GIEOgEPrcDjwG5wCjgERH5sZ3BlHJSxYE2vAbG5aY6HeUwM8b4irA2SUSGQHscrgCO8+uc+098Q9XusimXUodIjk8O6fn21LYAMD4MmyMyk+MpzknWK+EIEWhzxD78btIAEoGKoXYSkYUiUiIipSJySz/rE0XkKWv9+yIywVoeLyKPicinIrJFRG4NMKeKUC9f8TIvX/FyyM7XOzwtHJsjAGaNyWTzvganY6ggGLQIi8h/i8hvgAZgs4g8KiKPAJ8xxDhhEYkF7gcWATOBpSIys89m1wD1xpgpwK+Ae6zllwKJxphjgBOB63oLtFKhsKe2leT4WPLSE52O0q9ZYzPYXdtKU3uX01HUERqqOWKD9eeHwCq/5WsCOPZcoNQYsxNARFYCFwKf+21zIfAT6/UzwG9FRPBNGpQqInFAMtAJ6P+9othP//FTAG4747aQnG9PbQvjc1N8z3YLQ7MKfe3CW/Y3MXdijsNp1JEY6vFGj/W+FpEEYJr1tsQYM9Q/wYVAmd/7cuDkgbYxxnSLSAO+zr9n8BXo/UAK8G/GmLq+JxCR5cBygHHjxg0RR7nZG7veAEJZhFuZOCr8OuV6zRzju2V5874GLcIuF+joiAXAdnzNC78DtonI6fbFYi7QA4wFJgI/EJFJfTcyxjxojJljjJmTl5dnYxwVTbxew9661rDslOtVkJFIbmqCds5FgEBHR9wLnGOMKQEQkWnAn/G11w6kAij2e1/E4Z15vduUW00PmUAt8HXgFetqu0pE3gXmADsDzKvUiFU1ddDR7Q3L4Wm9RETvnIsQgY6OiO8twADGmG1A/BD7rAemishEqyljCfB8n22eB662Xi8G3jTGGGAv1gxtIpKKb9a2rQFmVeqIHByeFqYjI3rNGpvJ9qomOru9TkdRRyDQK+EPReQPfPFkjSv4otOuX1Yb7434HoMUCzxsjNksIncCG4wxz+O7/flxESkF6vAVavA1ezwiIpvxPdXjEWPMJ8P5xlRkyU0J3cO99/QOTwvj5gjwjZDo6jFs8zRxdKFOa+lWgRbh64HvAP9qvf8nvrbhQRljVgOr+yy73e91O77haH33a+5vuYpef73sryE7157aFmJjhLFZob1BZLiOsQrvJ+UNWoRdbMgibI33/dgYcxTwS/sjKRV6J59yKpUej/XmKsgpZsq0L2Zr9fSuc5Cn0sP4yVMPXXjBXdz6y//l1g0rGV1QwPtr33EmnBqxIYuwMabHuuttnDFmbyhCKdXXra/7bpq8+yt323L8So+HFQ+8BMCfP9hLUnwsF1nvAW6+4ARbzjscPV7vwYy9nttUQXPO6Vx57Te47/qvOpRMHYlAmyOy8d0x9wHQ0rvQGHOBLamU6mNd+bqQnauhrYuCjKShNwwDozOSeL+2TjvnXCzQIhyaEfJKOay9q4eObi9ZyUMN/gkPo61/LDyNOrewWw1ahEUkCV+n3BTgU+AhY4w+a1tFrANtvhtBM1PcUYQLMn1FuFKLsGsNNU74MXw3SXyKbyKee21PpJSDGlqtIuySK+Hk+Fgyk+P1StjFhmqOmGnNZIaIPAR8YH8kpQ5XlFEUkvMcaPM9MMYtRRh8TRLlB1qdjqFGaKgifHCSHuvmC5vjKNW/Jy5+YuiNgqC+tYv0pDjiYwN+8pfjRmcmUeJpgmQdK+xGQxXh40Sk9+Z0AZKt9wIYY0yGremUCrH6lk6yUxKcjjEsBRnWnMc5450NokZkqKksY0MVRKnBrHhlBQD3LbzPtnMYYzjQ2sVRY9wxPK1XXloiMQI9WoRdKdAhako5alPlJtvP0drZQ2eP13VXwnGxMeSlJ1KpRdiV3NPwpZTN6lt9nXLZLhme5m90RhLkFNPdozdtuI0WYaUs9dbwNLddCQOMyUxG4hL5fL/OL+w2WoSVshxo7SQ2RkhPcl8r3dgsXzv2+t31DidRw6VFWLnCtNxpTMudNvSGR6C+tYuslPiwfbjnYNKT4jEttWzYfdijGFWYc98/+SoqPXj+g7afo76lk1Fh+oj7gFTvZP3uMRhjXPkPSbTSK2GlACSWhvYuV3bKHVSzk5rmTnbX6t1zbqJFWLnC8heWs/yF5fadIDUXY9zZKXdQje85uOt3aZOEm2gRVq6wrXYb22q32XeC9DzA5UW4yUN2SjzrtV3YVbQIKwWQng9AlpubI4A5E3K0CLuMFmGlANLzSY6PJSne3XfqnzQhm921rVQ16dSWbqFFWCmA9Hx3d8pZTpqQA8AGHS/sGlqElSvMHj2b2aNn23eC9HyyU13cHmyZNTaTpPgYPtDOOdfQccLKFeycPa2hrQtJSnd9ezBAQlwMJ47P5n0twq6hV8Iq6m33NAGQm+riGzX8zJ+Uy5b9jdS1dDodRQVAi7ByhSufvZIrn73SlmNvrewtwu5vjgCYPzkXgPd31jqcRAVCi7ByhfLGcsoby2059jZPE6ar3ZUT9/Tn2KIsUhJiWbtDi7AbaBFWUW9rZRM07I+Y+RbiY2M4aUIO6/RK2BVsLcIislBESkSkVERu6Wd9oog8Za1/X0Qm+K07VkTWichmEflURNz1zBnlCsYYtnl8RTiSnDI5l9KqZh0v7AK2FWERiQXuBxYBM4GlIjKzz2bXAPXGmCnAr4B7rH3jgCeA640xs4AF+D35WalgqWrq4EBrFzRGVhHubRdep00SYc/OK+G5QKkxZqcxphNYCVzYZ5sLgces188AZ4nv/4TnAJ8YYz4GMMbUGmN6bMyqwtz8ovnML5of9OOWWJ1ykXYlPGtsJulJcbynTRJhz86eiEKgzO99OXDyQNsYY7pFpAHIBaYBRkReBfKAlcaYn/c9gYgsB5YDjBs3LujfgAofd3/lbluOG6lFODZGOHlirnbOuUC4dszFAacCV1h/XiQiZ/XdyBjzoDFmjjFmTl5eXqgzqghQ4mkiLz0ROlucjhJ08yfnsqe2lX0H2pyOogZhZxGuAIr93hdZy/rdxmoHzgRq8V01v22MqTHGtAKrgRNszKrC3CVPX8IlT18S9OOWVDYxvSA96McNB6dou7Ar2FmE1wNTRWSiiCQAS4Dn+2zzPHC19Xox8KYxxgCvAseISIpVnM8APrcxqwpzta211LYGt5j0eA3bq5qYFqFFeHpBOtkp8dokEeZsaxO22nhvxFdQY4GHjTGbReROYIMx5nngIeBxESkF6vAVaowx9SLyS3yF3ACrjTEv2ZVVRae9da20d3k5anRkFGFPpYfxk6ceunDeMp6pHsczty5mdEEB7699x5lwakC23iJkjFmNrynBf9ntfq/bgUsH2PcJfMPUlLJFb6fctAgpwj1eLyseOPRa5ePyA6wpqebqX67i0e9f5FAyNZhw7ZhTynbbrIl7phWkOZzEPsXZKQCU12vnXLiKjJvlVcQ7a+Jhg2OOWEllE+NyUkhJiNy/Btkp8aQmxFJWr09gDleR+9unIsptZ9wW9GNuqWxkeoQ0RQxERCjKSWFvrRbhcKXNESoqNXd0s6umhWMKM52OYrtxOSm0dfVA5hino6h+aBFWrrDoyUUsenJR0I73+b5GjIGjCzOCdsxwVZyd7HuRP93ZIKpf2hyhXKGtK7gdS59VNABwdBRcCacnxZOdEk9dwTSno6h+6JWwikqfVTSQn55Ifnp0zJA6LicF8ibT0a3zYIUbLcIqKn22ryEq2oN7jctJQeIS+WjvAaejqD60CKuo09rZTWlVM7OiqAgXZidjvD28W1rjdBTVhxZh5QrnTTuP86adF5Rjbd7XiNcQVVfCiXGxULeXf27XIhxutGNOucJNp9wUtGNtsv5LPrs4K2jHdAVPCZ+UT6ShrYvM5Hin0yiLXgmrqLOp7ABF2cm+eYSjSdU2vAZ92kaY0SKsXGHBowtY8OiCoBzro731HD8uOyjHcpXaPaQkxPKONkmEFS3CKqpUNbazr6E9+poiAEwP8yblaudcmNEirKLKR2UHgChsD7Z8acoodta0UKGPPAobWoRVVNm4p574WGHW2Mi/Xbk/p04ZBaBXw2FEi7CKKh/sruPYoiyS4mOdjuKIaQVp5KUnartwGNEhasoVLpt12REfo62zh0/LG7j2tElBSOROIsKpU0bx9rZqvF5DTIw4HSnqaRFWrnDDSTcc8TE+Kqun22uYOzEKR0b4+dKUUaz6qIISTxMzxkRns0w40eYI5QqtXa20dh3ZxOTrd9UjAieOzwlSKnfqbRfWJonwoEVYucK5T57LuU+ee0THWL+7jukF6VF/t9jozCSm5KfxjnbOhQVtjlBRoaO7hw176lhy0jinozjGU+lh/OSpvjezL2b7xJMZP/Uq8PqmtxxdUMD7a99xMGF00iKsosLGPQdo7/JyyuRcp6M4psfrZcUDLwGws7qZFz7ZzyX/byVF1hOZ77v+q07Gi1raHKGiwtodNcQIzIviIuyvMDsZEdhbpw8AdZoWYRUV3i2t4diiLDKSors9uFdiXCyjM5Ioq9M755ymzRHKFZbNXjbo+pNPOZVKj6f/lXGJyNfu5jtfnhr8YC5WnJPC+l11dHT1kBilN6+EAy3CyhWGKsKVHs/B9s6+dlQ38+In+zl1Sp4NydxrXE4KH+yqo6y+jSn5aU7HiVraHKFcoaa1hprWkQ2p2l3TgulqZ86E6L5Jo6/RGUnExwpl2i7sKFuLsIgsFJESESkVkVv6WZ8oIk9Z698XkQl91o8TkWYRCd5jFZQrLX56MYufXjzs/Ywx7K5tBc9W4mP1msNfbIxQmJXM3notwk6y7bdSRGKB+4FFwExgqYjM7LPZNUC9MWYK8Cvgnj7rfwm8bFdGFflqmjtp7uiG/VucjhKWxuWkcKC1i8b2LqejRC07Lw3mAqXGmJ3GmE5gJXBhn20uBB6zXj8DnCUiAiAiXwN2AZttzKgi3K7aFt+LSi3C/RmX4xsjrE0SzrGzCBcCZX7vy61l/W5jjOkGGoBcEUkDfgjcMdgJRGS5iGwQkQ3V1dVBC64ix87qZvLTE6G90ekoYSknNYGUhFgdL+ygcG0k+wnwK2NM82AbGWMeNMbMMcbMycvTnm91qKb2LjyNHdrzPwgRYVxOijVeWKe1dIKdQ9QqgGK/90XWsv62KReROCATqAVOBhaLyM+BLMArIu3GmN/amFeFsW/P+faw99lZ7WuKmJyXxrvBDhRBxuWksLWyCTLHOB0lKtlZhNcDU0VkIr5iuwT4ep9tngeuBtYBi4E3jTEGOK13AxH5CdCsBTi6XX705cPep7S6meyUeHJSE2xIFDmKrbkjKJjmbJAoZVtzhNXGeyPwKrAFeNoYs1lE7hSRC6zNHsLXBlwKfB84bBibUgBlDWWUNZQNvaGlpaObCr0JISBpSXHkpCRAwXSno0QlW++YM8asBlb3WXa73+t24NIhjvETW8IpV/nGqm8AsGbZmoC23+ZpwgDTC9LtCxVBxuWkUNs4mY7uHhLj9BbmUArXjjmljkiJp4m8tERy0xKdjuIKxbnJSFwCH+6udzpK1NEirCJOfWsnnsYOpo/Wq+BAFWWlYLzd/GObDvUMNS3CKuKUVDYBvse7q8AkxMVA9U4twg7QIqwiijGGksomirKSSde5g4encgtbK5vwNLY7nSSqaBFWrvCD+T/gB/N/MOR2nqYODrR1aVPESFRuBdCr4RDT+YSVK5w//fyAtivZ30SsiA5NG4nG/RRkJPKPbdVcNqd46O1VUOiVsHKFkpoSSmpKBt2mu8fLlspGJuelkqRPihiR06fm8c72Gnq8xukoUUOLsHKF6168jutevG7QbUqrm+no9jKrMDNEqSLPGdPzaGjr4uPyA05HiRpahFXE+KyikczkeIqzk52O4lqnThlFjMA/SrRdOFS0CKvIkJZHxYE2Zo3NwJqSWo1AVkoCs4uzWKOdcyGjHXMqMkyaT4zAzDEZ/a72VHoYP3ngpy17BnpScxQ6c3o+v3x9G1VN7eSnJzkdJ+JpEVau19nthfEnMXFUKqmJ/f9K93i9Az6NGeDmC06wK57rfGVmAff+fRtvba3i8pPGOR0n4mkRVq7w49N/POC6VzdXIknpHK0dckFx1Oh0CrOSeX2LFuFQ0CKsXOErk74y4Lo/rtuNaa5mfM6UECaKXCLCWTPy+cuGctq7enS4n820Y065wqbKTWyq3HTY8s37Gli/ux5K39EOuSA6a0YBbV09rNtR63SUiKdFWLnCildWsOKVFYct/+PaPSTHx8LuD0IfKoLNm5RDakIsf9+iHZZ20yKsXKu+pZPnNlVw0QmF0NXmdJyIkhgXy+nT8nhzSxW+J44pu2gRVq719IYyOrq9XDV/vNNRItJZMwqobGxn875Gp6NENC3CypV6vIbH39vDvEk5HDW6/7HB6sicOT2PGIHXPtcmCTtpEVau9ObWKsrr27h6/gSno0SM3htaer9OOO5oejzb+fWqfzJ+8lROPuVUpyNGJB2iplzhZ2f97JD3j63dzZjMJM6eWeBQosjT3w0tn5Qf4K2Saq74xbM8cdPFDiWLbHolrFzhlOJTOKX4FABKq5p4p7SGK+eNJy5Wf4XtNCU/DQG2e5qdjhKx9DdYucLasrWsLVsLwEPv7CYhLobLT9KJx+2WkhBHUXYy2zxNTkeJWNocoVzhR2/8CIC/Ln6VZzeWc8kJhYzSx9mHxLSCdN7YWgWZhU5HiUh6Jaxc5Yn39tLR7eWaUyc5HSVqTM5LQwQoPt7pKBFJi7ByDa8xPP7ebr58VL4+Qy6EkhNiKc5OgeLZeuOGDbQIK9eoae6kprmTa0+d6HSUqDOtIA1JG8WmsgNOR4k4WoSVa+xvaGPmmAzmT851OkrUmZKfhunu4OkNZU5HiTi2FmERWSgiJSJSKiK39LM+UUSesta/LyITrOVni8iHIvKp9eeX7cypwt/+l9JJbrmGzc/dz4Qp0w65qWD85Kn6ZAybJcbFQtkmnt+0j+aObqfjRBTbRkeISCxwP3A2UA6sF5HnjTGf+212DVBvjJkiIkuAe4DLgRrgfGPMPhE5GngV0K7ZKNae9mVyEqax7JZziI05fMpKfTJGCOxaR8vEk3nh430snauTvQeLnVfCc4FSY8xOY0wnsBK4sM82FwKPWa+fAc4SETHGfGSM2Wct3wwki4iOR4pSm/c10D6mlYK8kn4LsAqR2t1MzU9j5Qd7nU4SUewswoWAfwNSOYdfzR7cxhjTDTQAfRv8LgE2GmM6+p5ARJaLyAYR2VBdrU+HjVS/e2sHDbF/YlfnE05HiXpL5o7j4/IGPteZ1YImrDvmRGQWviaK6/pbb4x50BgzxxgzJy8vL7ThVEiUVjWx+rP90NZAjD45w3EXH19IQlwMK9fr1XCw2FmEKwD/+0qLrGX9biMicUAmUGu9LwJWAVcZY3bYmFOFsd++Wep7ckbbAaejKCA7NYFFR49m1cYKmtq7nI4TEewswuuBqSIyUUQSgCXA8322eR642nq9GHjTGGNEJAt4CbjFGPOujRlVGNtV08LzH+/jynnjwdvjdBxl+daXJtLU0c3KD3S4WjDYVoStNt4b8Y1s2AI8bYzZLCJ3isgF1mYPAbkiUgp8H+gdxnYjMAW4XUQ2WV/5dmVV4el3b5USHxvDtafpzRnh5LjiLE6emMPD7+6iq8frdBzXs3UCH2PMamB1n2W3+71uBy7tZ7+7gLvszKbC266aFlZ9VMGV88aTn55Eznv5XLroDqdjKcv1Z0zmm4+uZ9XGCi7T2eyOSFh3zKnodc/LW0mMi+GGMycDEN+YQH6yTtoTLhZMz+PYokx+/cZ2Orv1avhIaBFWYWf97jpe2VzJ9WdMJj89CYDWohY2173pcDLVS0T4/tnTqDjQprcyHyEtwiqsGGO466UtFGQkcu1pX1z5Ns2sZ82+RxxMpvo6Y1oecyfkcN/r22jUkRIjpkVYhZUXP9nPx2UH+ME500lOiHU6jhqEiHDbeTOpbenkt2+WOh3HtfTJGipstHR0c/fqLRw1Op1LTihyOo7qo/dpzIc58XJ+39PDM7/6dzb+/dnQB3M5LcIqbNz72jb2N7bz318/XueICEP9PY0ZoK2rhyfe20PtpH+hq8dLvD58dVj001JhYf3uOh5du4srTx7PieNznI6jhiE5PpYzp+cj2UXc+9o2p+O4jhZh5bjG9i5WrNxEUXYKP1x0VL/b5L5TwNen/jzEyVSgpuSnYXas5YF/7OD1z3Vu5+HQIqwc5fUa/u9fPqaysZ1fXT6btMT+W8jiWuPJThwT4nRqWDY9y6yxGXxv5Ud8VtHgdBrX0CKsHPWbN7fz6mYP3Ruf5eIvzz3siRm9X/uyyvmoZvXQB1TO8Xbz8LKTyEyOZ9kj6ymtanI6kStox5xyzBPv7eG+17djdn/A977/fxG5ecBtv79yBmsr/8zxo84NYUI1XAUZSfzxmrks/d/3ufz37/HYt+ZydGGm07HCml4JK0c8+u4ubvvbZ5x1VD5sWInoXMERY0p+Ok8tn0dSfCyLH1jL3zZVYIxxOlbY0iKsQqqz28tPnt/MT174nLNnFHD/FSeA0bkHIs2kvDSe+86XOHpsJt9buYkbntxITfNhD8dRaBFWIfRZRQMX/8+7PLp2N9/80gR+d8UJJMXrXXGRKi89kZXL5/HDhUfxxpYqzvnV2/zp/b06/WUf2iasbLfN08Tv3irlbx/vIzc1kQeuPIGFR+tIh2gQFxvDtxdM5qwZ+fzo2U/50apP+f3bO1jxlamcd+xYvbEDLcLKJjuqm3n9cw+vfe7hwz31JMfHsvy0Sdxw5hQyk+OHfbzYVUlcfcmvbUiqQmFaQTp/uX4+b5VU8V+vbuPfnvqYn79SwlXzJ7B0bjFZKQlOR3SMFmEVFN09XjbuPcDrWzy8/rmHnTUtAMwck8Eti47isjnF5KSO/C+atAlp8XonnZuJCF8+qoAF0/J5c2sVD7+7i3te2cpv3tjOxScUcuW88cwYk+F0zJDTIqxGrLPby7s7avj23X+gLXsykpiG8XZD1XbYtxn2fcbmtgP8s7aWu3NzBzyOxzP0HVbeY7r4oOpZ5uZfHMxvQQXRgBP8WOpqa8np+3uQOYbWqQt4ov0Ennx/L12eUuL2roeyj6Cn87BjjC4o4P217wQ7uqO0CKthMcawbmctz31UwaubPTS0dWFypjFjXD6TRqUyLjeFxLgZwAUH97n5ghP6nfjFf/1QvMd0s75qlRbhMDbQBD+9Bvs9aOvqYev+Rt6oq0BOWkrCvCs4anQ6M8ZmUJCeeHAI433Xf9WW7E7SIqwC0tDWxV8/LOfJ9/ewo7qF9MQ4zp5ZwFePHcO3LjiDf/ld3wdpKxW45PhYjh+XzZ9v/D7fe+xdPtvXwOb9jXxS0UBmcjzTC9KZVpDmdExbaBFWg9ruaeLhd3ez6qNy2ru8zC7O4t5Lj+Orx475YniZPo5eBVFhdjKF2cmcMa2H0upmtlU2sX53HR/sroOFP+L2v33GqVNGMX9yLulJw+/kDTdahNVhjDG8vb2Gh97ZxdvbqkmMi+Frswv5xvzxeguqCpmk+FiOHpvJ0WMzaenoZntVM2v+sZm/bCjnj+v2EBsjzBiTzqwxmRxdmMGswkymF6STOsAkUOHKXWltcPIpp1I5SMdQJHYE9NX7Gazctw+S0lhyw0NI5hhMWwPseIf9G15kZWoiKwfYP5CONaWORGpiHLOLs3jqmbsoGDMWcifQXTCdT/eN49PsIiTxi6aKgoxEJo5KZeKoNCaOSrH+TGVcTgqnnX76oH/f++089GNHPYj6Ilzp8QzamRCJHQF9VbbCnB88Ruw91+M1hvzCCRw/LoupBZOJi5nDzRf8kTse3zjg/oF0rB2p2KeT+D+XPmj7eVR46/F6WfE/h/Y/GGNo6uimuqmDF578X05bdj27alp4dXMldS1fjLCIEeg5/lomTZxATkoC2anx5KQmkJOSQKLVtDZUJ7Id9SDqi3C0qmvp5KVP9/PMh+XIwlv5qKyehLgY0pLiWDq3OOwm1JFuISE22ekYKgyJCBlJ8WQkxcPW1/nFpf9zcF1Daxe7alvYVdPMzuoWfvPYRzS2F7K3tpUev0mFUhJiyUtLJGP+5WyvaqIgPYn0pLiQ/D3QIhwlunq8bN7XyLodtbyxxcPGvfV4DRw1Oh2zaRXXfPf75K5OBAi7AgzgPb6Ld/c/yZfGXOF0FOUimSnxzE7JYnZxFgC/+fZjXLnsMrzG0NjWRV1LJ3WtndS1dFLd1EHanAtZ/WklAEnxMRRkJDE2K5nCrGQKMhJtyahF2AWG025tjKGhrYsd1c2UVDazzdPEk6v/QWfaaCTO90tk6susmyk+ZcuBCjweDyk/GHgu33DgndHNptpXtAirQQ11w0hv/0WMCFkpCWSlJDDJb/3NF53M9x/9B56mdqoaO6hsbGfdjloA38Nnz7gx6A8zjeoi/MNnPoHTruPpDWV09XiJESE2RoiLFVIS4khNiIUpp/H2tmom5aUyNjOZGAeeAtzbbm2MobWzh+aObt9Xu+/P9W+/xtIH36OysZ39DW20d30xS1VqQiydPYbZE/IP/ouemjgV+PLBbULRpqtUKARyw8jgB+hidGYSozOTDi5q6+xhX0MbFQfa2Li7IeiTDtlahEVkIfBrIBb4gzHmP/usTwT+CJwI1AKXG2N2W+tuBa4BeoB/Nca8Gux8+xraID6FuBghOT4erzF4je923P2tbbR09iDHX8JVD38A+P57MnFUGpPzUpmU5/tzcl4ak/JSSUk4so+yu8dLTXMnnsZ29jf4imllg+81C77LI+/uoqWj55B2LPB1NjBqMl09XmaOzeCso/IZnZnExFGpTCtIpzArmYlTp7HgsoF/MZVSA0tOiGVyXhqT89LYeO/jwO1BPb5tRVhEYoH7gbOBcmC9iDxvjPncb7NrgHpjzBQRWQLcA1wuIjOBJcAsYCzwuohMM8YE9a6Ax685mfE/u5KLByhQxhh+/W9f56kXXmNndQs7q5vZWdPCpxUNrP50P16/ejjW+tczNy2RUWkJ5KQmkBwfS2xMDHExvivsjm4vrZ2+q9eWjm7qWjrxNHbgaWynprnjkOMBJMTFMCYzCYxhTFYyaYlxpCfGkZYUR1qi7yslIZZff/s8nvnv7cH8aJRSIWLnlfBcoNQYsxNARFYCFwL+RfhC4CfW62eA34qvV+hCYKUxpgPYJSKl1vHW2Zj3MCIC7Y3Mm5TLvEmHjh3s6O5hT20rO6qa2VHt63mtauqgrK6VTWUHqGvppKdvVQVEIC0hjpTEWHJSEynISGTmmAwKMhLJz0iiICOJMZm+r5zUBESE8ZOXs3CpXskqFYnErmc/ichiYKEx5lrr/TeAk40xN/pt85m1Tbn1fgdwMr7C/J4x5glr+UPAy8aYZ/qcYzmw3Ho7HSgJ4rcwCqgJ4vGCQTMFLhxzhWMmCM9c4ZgJRp6rxhizsL8Vru6YM8Y8CNgygl9ENhhj5thx7JHSTIELx1zhmAnCM1c4ZgJ7ctn5bJEKoNjvfZG1rN9tRCQOyMTXQRfIvkop5Xp2FuH1wFQRmSgiCfg62vrOd/g8cLX1ejHwpvG1jzwPLBGRRBGZCEwFPrAxq1JKOcK25ghjTLeI3Ai8im+I2sPGmM0iciewwRjzPPAQ8LjV8VaHr1Bjbfc0vk68buA7wR4ZEYBwnKhAMwUuHHOFYyYIz1zhmAlsyGVbx5xSSqmh6fOmlVLKQVqElVLKQVFVhEUkR0T+LiLbrT+zB9juFRE5ICIv9lk+UUTeF5FSEXnK6nDE6kB8ylr+vohMsCnX1dY220XkamtZuohs8vuqEZH7rHXLRKTab921ochkLV8jIiV+5863ljv5WaWIyEsislVENovIf/ptP+zPSkQWWt9jqYjc0s/6Ab9XEbnVWl4iIv8S6DHtyiQiZ4vIhyLyqfXnl/326fdnGaJcE0Skze/cD/jtc6KVt1REfiMyvOn/jiDTFX3+znlFZLa1bviflTEmar6AnwO3WK9vAe4ZYLuzgPOBF/ssfxpYYr1+APi29foG4AHr9RLgqWDnAnKAndaf2dbr7H62+xA43Xq9DPitXZ/VYJmANcCcfvZx7LMCUoAzrW0SgH8Ci0byWeHrbN4BTLKO9TEwM5DvFZhpbZ8ITLSOExvIMW3MdDww1np9NFDht0+/P8sQ5ZoAfDbAcT8A5gECvNz7s7Q7U59tjgF2HMlnNaIP1a1f+O6oG2O9HgOUDLLtAvyKsPWDrgHirPfzgVet168C863XcdZ2EsxcwFLg937vfw8s7bPNNKCs99wcWRE+okwD/TKGy2dlLf818H9G8ln5//yt97cCtwbyvfbdtne7QI5pV6Y+2wi+0UqJg/0sQ/RZTaCfImz97LcO9DMP4Wf1M+D/+b0f9mcVVc0RQIExZr/1uhIoGMa+ucABY0y39b4cKLReF+IrfljrG6ztg5nr4Dn6OX+v3n+t/Ye8XCIin4jIMyJSTOCCkekR679kt/n9VzEsPisRycL3v503/BYP57MK5Ocx0Pc60L6BHNOuTP4uATYa39wtvfr7WYYq10QR+UhE/iEip/ltXz7EMe3M1Oty4M99lg3rs3L1bcv9EZHXgdH9rPp3/zfGGCMiIRufF6JcS4Bv+L1/AfizMaZDRK4DHsNvImGbM11hjKkQkXTgr1auPwayo92flfjuzvwz8BtjTTDFEJ9VtBCRWfhmMzzHb/GIf5ZBsB8YZ4ypFZETgeesjI4TkZOBVmPMZ36Lh/1ZRVwRNsZ8ZaB1IuIRkTHGmP0iMgaoGsaha4EsEYmz/lX0v5W69zbrcjn09utg5qrA10TSqwjff316j3EcvqaSD/3O6Z/hD/jaU0OSyRhTYf3ZJCJ/wjcL3h8Jg88K34D77caY+/zOOehnNcA5Ar0tv+/3Oti+R3K7/pFkQkSKgFXAVcaYHb07DPKztD2X9b+6Duv8H4pvkq9p1vZFQxzTlkx+65fQ5yp4JJ9VtDVH+N8mfTXwt0B3tH4Z3sJ3e3Xf/Qe6/TqYuV4FzhGRbPGNCDjHWtZrKX1+Iawi1esCYEsoMolInIiMsjLEA+cBvVcLjn5WInIXvr9MK/x3GMFnZcdt+YEc05ZMVvPMS/g6Pd/t3XiIn2UocuWJb25yRGQSvs9qp9Uk1Sgi86z/8l/FMP4+H0kmK0sMcBmwsnfjEX9WI21sd+MXvvacN4DtwOtAjrV8Dr4nf/Ru90+gGmjD11b0L9bySfj+spQCf+GLjosk632ptX6STbm+ZZ2jFPhmn2PsBI7qs+xuYDO+nt+3+q63KxOQim+UxifW+X8NxDr9WeG72jH4Cuwm6+vakX5WwLnANny97P9uLbsTuGCo7xVf08oOfB2NiwY75jA/nxFlAn4MtPh9LpuA/MF+liHKdYl13k3ARuB8v2POwVfkdgC/ZRgdvEH4+S3AN92u//FG9FnpbctKKeWgaGuOUEqpsKJFWCmlHKRFWCmlHKRFWCmlHKRFWCmlHKRFWEUlEemxbi39TET+IiIpQ2y/RkTC7sGTyv20CKto1WaMmW2MORroBK53OpCKTlqElfLdnDNFRBZYV7zPiG/O4SdHMFmNUsOiRVhFNWtOgEXAp9ai4/Hd0jwT3x2SX3ImmYoWWoRVtEoWkU3ABmAvvid/A3xgjCk3xnjx3So7wZF0KmpE3CxqSgWozRgz23+B1fLgP4duD/p3RNlMr4SVUspBWoSVCtxLIlJuff3F6TAqMugsakop5SC9ElZKKQdpEVZKKQdpEVZKKQdpEVZKKQdpEVZKKQdpEVZKKQdpEVZKKQf9f8bV/uarzngGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = sns.displot(PnL, kde=True, stat=\"probability\").axes.flatten()\n",
    "fig[0].axvline(var_1p, c=\"r\", linestyle=\"-\", label=\"1% VaR\")\n",
    "fig[0].axvline(var_5p, c=\"g\", linestyle=\"--\", label=\"5% VaR\")\n",
    "fig[0].legend()\n",
    "fig[0].set(xlabel='PnL', ylabel='Probability')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "40a6b933-737f-48af-9749-9a25a13157a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2646a37781a843aa9ed644dbc7c05b63",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/5000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "735ac97e78754e8da747b4e80a575b8e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/5000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bootstrap 99% VaR: -0.06626183287860131\n",
      "Confidence intervals: [-0.09966348 -0.05024794]\n",
      "Bootstrap 95% VaR: -0.035440180885582526\n",
      "Confidence intervals: [-0.0478463  -0.02684057]\n"
     ]
    }
   ],
   "source": [
    "# Bootstrap VaR\n",
    "def VaR1(data):\n",
    "    return np.percentile(data, 1, interpolation=\"lower\")\n",
    "\n",
    "\n",
    "def VaR5(data):\n",
    "    return np.percentile(data, 5, interpolation=\"lower\")\n",
    "\n",
    "\n",
    "def bootstrap(data, func):\n",
    "    sample = np.random.choice(data, len(data))\n",
    "    return func(sample)\n",
    "\n",
    "\n",
    "def draw_bs_reps(data, func, size):\n",
    "    bs_reps = np.empty(size)\n",
    "    for i in tqdm(range(size)):\n",
    "        bs_reps[i] = bootstrap(data, func)\n",
    "    return bs_reps\n",
    "\n",
    "\n",
    "bs_var_1 = draw_bs_reps(PnL, VaR1, 5000)\n",
    "bs_var_5 = draw_bs_reps(PnL, VaR5, 5000)\n",
    "\n",
    "print(f\"Bootstrap 99% VaR: {np.mean(bs_var_1)}\")\n",
    "print(f\"Confidence intervals: {np.percentile(bs_var_1, [2.5, 97.5])}\")\n",
    "\n",
    "print(f\"Bootstrap 95% VaR: {np.mean(bs_var_5)}\")\n",
    "print(f\"Confidence intervals: {np.percentile(bs_var_5, [2.5, 97.5])}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "id": "55c96f8b-0f86-4e29-8c96-aa554fd34f78",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.08570424750973445"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "-0.0522112397056073"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Conditional VaR - Expected tail loss - Expected shortfall\n",
    "exp_sf_1p = PnL[PnL < var_1p].mean()\n",
    "exp_sf_5p = PnL[PnL < var_5p].mean()\n",
    "\n",
    "exp_sf_1p\n",
    "exp_sf_5p"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}