{ "cells": [ { "cell_type": "markdown", "id": "6f883584-ca8b-4193-82bd-eca78426165e", "metadata": {}, "source": [ "# Méthode de Monte Carlo\n", "\n", "Tout d'abord, on illustre numériquement les deux résultats probabilistes sur lesquels reposent la méthode dite de Monte Carlo: \n", "\n", "- la loi forte de grands nombres,\n", "- le théorème central limit (TCL).\n", "\n", "On considère ensuite un premier exemple d'estimateur de Monte Carlo et l'importance de l'intervalle de confiance (IC) dans lequel se trouve la valeur recherchée avec probabilité grande (0.95). \n", "\n", "Enfin on applique la méthode de Monte Carlo à un exemple multidimensionnel où on illustre l'efficacité de 2 méthodes de réduction de variance: \n", "\n", "- variables antithétiques, \n", "- variable de contrôle. " ] }, { "cell_type": "code", "execution_count": 1, "id": "2e046f21-4713-4bfc-b2d0-32073d10d250", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "sns.set_theme() \n", "from numpy.random import default_rng\n", "rng = default_rng()" ] }, { "cell_type": "markdown", "id": "dc0d798f-524f-46de-aa81-2e7107a546d5", "metadata": {}, "source": [ "## Illustration de la loi des grands nombres\n", "\n", "Soit $(X_n)_{n \\ge 1}$ une suite de variables aléatoires _i.i.d._ de carré intégrable. On définit les suites $(m_n)_{n \\ge 1}$ et $(\\sigma_n^2)_{n \\ge 2}$ (non définie pour $n = 1$) de la façon suivante\n", "\n", "$$\n", " m_n = \\frac{1}{n} \\sum_{k=1}^n X_k \\qquad \\text{et} \\qquad \n", " \\sigma_n^2 = \\frac{1}{n-1} \\sum_{k=1}^n (X_k - m_n)^2 \\quad \\text{pour} \\; n \\ge 2\n", "$$\n", "et on veut illustrer la Loi Forte des Grands Nombres et le Théorème Central Limite (étendu en utilisant le lemme de Slutsky pour remplacer $\\sigma^2 = \\mathrm{var}(X_1)$ par l'estimateur $\\sigma_n^2$) c'est à dire les convergences\n", "\n", "$$\n", " m_n \\xrightarrow{p.s.} m \\qquad \\text{et} \\qquad \n", " \\sqrt{n} \\Bigl(\\frac{m_n - m}{\\sigma_n}\\Bigr) \\xrightarrow{\\mathcal{L}} \\mathcal{N}(0, 1).\n", "$$\n", "\n", "Plus précisément on construit l'intervalle de confiance (asymptotique) à 95% à partir du TCL c'est à dire\n", "\n", "$$\n", " \\text{pour $n$ grand} \\quad \\mathbf{P} \\biggl( m \\in \n", " \\biggl[\n", " m_n - \\frac{1.96 \\sigma_n}{\\sqrt{n}}, \n", " m_n + \\frac{1.96 \\sigma_n}{\\sqrt{n}}\n", " \\biggr] \\biggr) \\simeq 0.95\n", "$$" ] }, { "cell_type": "markdown", "id": "853c2b0a-e6d8-40f1-a3af-f1732ef4ea2d", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: LFGN loi uniforme\n", "\n", "Reproduire le tracé suivant où les points (les croix 'x') sont les réalisations $X_n$ (en fonction de $n$) de loi uniforme sur $[-4,8]$. La ligne bleue (couleur 'C0', première couleur de la palette utilisée) correspond à la moyenne $m$, la courbe orangée (couleur 'C1') correspond à la suite $m_n$ et les lignes grises correspondent aux bornes de l'intervalle de confiance. La zone de confiance en jaune s'obtient par la méthode `fill_between` de `ax`.\n", "\n", "![](img/tcl_unif.png)" ] }, { "cell_type": "code", "execution_count": 17, "id": "51db966d-5681-4222-8217-5f9e83d8a681", "metadata": { "tags": [ "aremplir" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\lione\\AppData\\Local\\Temp\\ipykernel_62128\\3030083443.py:10: RuntimeWarning: divide by zero encountered in divide\n", " s_n = np.sqrt(np.cumsum((uniforms - m) ** 2) / np.arange(0, n))\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAF/CAYAAADZxC9bAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAhIdJREFUeJztvQeUZGd1rr0rds558kxPDprRJEmghISFLmAjy7b8gzC+JhgvX6ML2IBZ4oJkLgvZxkpEX4MuP8nINvwyMrKFZYGQjcKMNKPJ0xN7pmc6TedY8fzr3ae/6qrqqu7qngrd1e+zVmlaVafO+eqcr07t99vJYVmWJYQQQgghhBCSAs5UNiKEEEIIIYQQQAFBCCGEEEIISRkKCEIIIYQQQkjKUEAQQgghhBBCUoYCghBCCCGEEJIyFBCEEEIIIYSQlKGAIIQQQgghhKQMBQQhhBBCCCEkZdySh6A3Xjic2/54Tqcj52Mg8wPOBWLgXCDRcD4QA+cCmQ9zAcd2OByLV0DgxPf2juTs+G63U6qqSmRwcFSCwXDOxkFyD+cCMXAukGg4H4iBc4HMl7lQXV0iLldqAoIhTIQQQgghhJCUoYAghBBCCCGEpAwFBCGEEEIIISRlKCAIIYQQQgghKUMBQQghhBBCCEkZCghCCCGEEEJIylBAEEIIIYQQQlKGAoIQQgghhBCSMhQQhBBCCCGEkJShgCCEEEIIIYSkDAUEIYQQQgghJGUoIAghhBBCCCEpQwFBCCGEEEIISRkKiBzT3n5Jzp49lethEEIIIYQQkhIUEDnm4sVW6ei4nOthEEIIIYQQkhIUEDnGsizx+wMSDAZyPRRCCCGEEEJmhAIixzidThUPo6OjuR4KIYQQQgghM0IBkWNcLpcEAgEZHx/L9VAIIYQQQgiZEQqIHBMMBsWywhIKhXI9FEIIIYQQQmaEAiLHIHwJeRDhcDjXQyGEEEIIIWRGKCByCISD7YGggCCEEEIIIQsDCogcgrClcNiigCCEEEIIIQsGCogcCwjkP1BAEEIIIYSQhQIFRA4JhYIR4UABQQghhBBCFgIUEDn3QFjicDgpIAghhBBCyIKAAiLnORBhbSYHbwQhhBBCCCHzHQqIHIIKTBAQLhcEBPtAEEIIIYSQ+Q8FRA55raVHnj5aJCe7CiQcpoAghBBCCCHzHwqIHPLkLy/LWMApb1z2iM/PECZCCCGEEDL/oYCYJ1wZDOR6CIQQQgghhMwIBcQ8oXOAHghCCCGEEDL/oYDIESjfGk3PEMu4EkIIIYSQ+Q8FRI7wBWKTpsf8FBCEEEIIIWT+QwGRI8Z8sQIiGLKmeCUIIYQQQgiZb1BA5IjxuKpLwbCwlCshhBBCCJn3UEDkiNHxWAERUgHBMCZCCCGEEDK/oYDIEaPjgSkhTOxGTQghhBBC5jsUEDli1Oef4oGggCCEEEIIIfMdCogcMToW54EIixw9ekh6erpzNiZCCCGEEEJmggIiR4z6Agk8EEHp7GzP2ZgIIYQQQgiZCQqIeeKBCFkOCQaDwkquhBBCCCFkPkMBkSNG4pOowyKBQIC9IAghhBBCyLyGAmIelXENBkMUEIQQQgghZHELiL/927+V3/u935t2m76+PvnTP/1T2bNnj+zdu1cefPBBGRsbk3xm1BcrIEQcEgyFNIwp5X2Mjkpvb0/ax0YIIYQQQkhOBMT3v/99efTRR2fc7r777pPW1lb59re/LY899pi88MIL8sADD0g+M+qbWrI1GApLMBgb2jQdR4++IS0tx9M8MkIIIYQQQpLjlgzQ2dkpn/vc5+SVV16RVatWTbvtgQMH5NVXX5VnnnlGmpub9bm/+Iu/kA9+8IPy8Y9/XBoaGmSxCIhA0BK/36dhTA6HY8Z9YLtAwJ/y9oQQQgghhMxLD8TRo0fF4/HIT37yE9m+ffu02+7fv1/q6uoi4gEgjAkG8WuvvSb5yngiARGyNIQpHA6ntA9sh7yJ2XgtCCGEEEIImXceiNtuu00fqXormpqaYp7zer1SWVkp7e1z74ngducuP9zlcsb8m4ixwFSREAxb6k0YGRmUysoqcTqn/wyWFdZHOBwUt7swDSMnuZgLZHHAuUCi4XwgBs4FshDnQkYExGxAsjQEQzwFBQXi8/nmtE+n0yFVVSWSa8rLixI+j1wHhCtN3b5SPB5LWltPi8jKGK9MPBAaHo9TAgGHFBW558XnJbOfC2TxwblAouF8IAbOBbKQ5kLOBURhYaH4/f4pz0M8FBcXz2mf4bAlg4OjkiugHHHxBwfHJIT6rHEMjkz9vCBsuWRsbFQcjlHp7OyR6urGpMdAz4jR0XEZG/PJlSsD4nTSAzEfmWkukMUD5wKJhvOBGDgXyHyZCzh2qt6PnAuIxsZGee6552Keg6Do7++X+vr6Oe83iM5sOQYXP9E4hpIIiGDYIaFQSMVTsvcaxsbGJRAIas4E/o7ftqurQ8OgvN6CNHwScrXMdD3J4oFzgUTD+UAMnAtkIc2FnAdZofdDR0eHlnE1oCoT2LVrl+QjY/7EvR6QFoHEaHgXZkqktkVGSJPNsX00eO/Zs6f1QQghhBBCyIIWEDB6u7u7ZXx8XP8fVZp27twpH/vYx+TQoUPy8ssvy2c/+1m566678raE61hUBSa3czIXor3fUuM/lEJDufHxMQmHQ5poHS8g4MFBeVfshxBCCCGEkAUtIFBZ6cYbb9S+DwAr6F/5yldk2bJl8vu///vy0Y9+VG6++ea8biQ3HtWFutg7KSCOtIWktdelwsDnG59+H+N2p26HwzllWwiKeFFBCCGEEEJIOsh4DsRDDz0U8/8QCidPnox5rqamRh5//HFZLIz7Jz0DRR6RwSj7f/+lElle2acemukaxI2MjKh4gAdiqoCAB4ICghBCCCGE5GEOxGIkOgeixJs4eSYUshOkk+5jbFRcLpc+4svdmvClmcKgCCGEEEIImS0UEDlgLCqEqSiBgCgrK5voMp3ci4A8B3gf7ByI2HwHvJZKGBQhhBBCCCGzhQIiB4yOTwqDkgRVVouKilQQJAtDMuLCFhCuSOWmaAGB9/v9Pg2DIoQQQgghJF1QQOTYA1FaOPUSGFGQzAOB5/tGwvJfZ5xysc+poUrR3oaRkaGJ7aYPgyKEEEIIIWS2UEDkWECUF0/NY0fvEMuCVyGx8Q9R8EKLU850O+Q/joc1KRs5EQAeh8HBQSkoKJgQIRQQhBBCCCEkfVBA5IDxqCTq8uKpSRB2hJNjGg9EUAbGJy8d9MjYmF3WFdWbELrk9RbOmEdBCCGEEELIbKGAyHEjubIEZZhMikSyHIh4r4Iv5IqELY2Ojuj74IGYLo+CEEIIIYSQuUABkVMPhCXFhVMFxC9PBCQQspJ6D8Z9/pj/D4RdMjw8rH9DSCD8ye1267/0QBBCCCGEkHRCAZED/IGw/utyiHg8UwXElWFLjnYUaDhSIvqHY5/3h9ALYlw9DgMD/dpgzm5A56AHghBCCCGEpBUKiBwQCE0ICKclbrdL7rl56ZRtWro9UxrEGQaGY5/3BZ0qHrD9wMCAeDyeyGtMoiaEEEIIIemEAiKXHginiMvllt3rq+XG1VO9DckawQ0Mx4Yw+UJIuA5q/gP+HQ+65FJfSMJhu+QrIYQQQggh6WJqDVGSNQ+EWwWESxvCeRNcCdNhGiVae3t7ZMWKVfp8f5yAGAs4tOLS8PCQ+AIhefqQR3zBgOxa5pZVSbwYhBBCCCGEzAV6IHJAIGhFQpiMgChIICBMh+nTp0/K5cttka7SVwZjRcGY335+ZGRYLvRaWtYVvNZWoCVdCSGEEEIISRcUEFkmHLYkFJ4QEA7bA4FHoQdJz7GgQRxCkKAb4IkwAqJ3KDavYdSPBnK2gDDixJAsj4IQQgghhJC5QAGRZfzByR4QbpeI02kLiIIEAmLUZ3sgUO41FAqqiAB9I5P7MB4Ip9MOY4ravX08vz8iPAghhBBCCLlamAORJcbHx7RbtMtbGnnOTqK2Q5jccEfEAc+CSYKGeICIQInWoXEIgsntAyERy+HSbcaDsfsx73O7JyszEUIIIYQQMlfogcgSR468IadOnYj1QEwICCMi4kEuA6oqQQSEwxACIRkYHpOwNVVsoJkcvBVIqI4G72EpV0IIIYQQki4oILIEwovQGdqUcAXwOtjiwaWN33YvH5/S3wGiAO/V3IlQSC5c6ki4f1/IJUNDA+IPxzqV/EETBkUIIYQQQsjVQwGRRZCLEAhGCQinnQMB7wNCk9Y3OOT6FcOR131Bh4YfGfEBAdHW3pVw3yHLI1VVNTIWiL2kY/4wPRCEEEIIISRtUEBkEXgRYpOobQ8EvA94aDUmd2wIk+2BCEU8EMNjk96EiuLJcKXBcRGX2yODY7EJ06M+S/r6ejP90QghhBBCyCKBAiLLHojYECY7BwKgihK8EIXeye19Iaf2cbBLuMKTEJCR8UlvQmP55LZ9I5YMjFoyUSE2Qtjhlba2C5EKToQQQgghhFwNFBBZFhAjo5N5Dh4XQpdsLwLEA0KZCt2TXgV/yKFdqEMhOwcCYmJ8omkcKHYMRf7uGwlL78jUcq1By63CA+VcCSGEEEIIuVooILIsIM61tkb+3+uZPP0QEhAQHpclTsdkDoTP59cQJngg0BRuPDApEgrdQSn22B4NiIfekUnvhsEfcqr3IRBIj4AYHR2V0dERLUvr88UmfRNCCCGEkPyHfSAyBHIXBgcHpLi4LC6EaTIEyYsYpgkQwmRZtheiwIPkZ1tA2PkPtjCAwR4VwaT5EuWFlowG7F4QF3umCghsDwGRLg/E4cMHtLEdetNB9DQ1LZHly1dFPCmEEEIIISS/oQciQxw6dEiOHTsS8xy8CJZjUrMVeCcFhAlhgpDwOicEQ9AhwaDd/8HlcsvIyIj4oiqyFnodUlIw+f/dQ1NDmHwTfSEQ/pSWKlKBgAQCQfWGIDTqwoXWtHk3CCGEEELI/IcCIkOMj48nNMDDUR2koz0Q1dU1kYZyHqed8BwKi4z77UZyXq9XjXdUZjIUepwxAsJQVzZ5jLGJkKd0eCAgQiAWUBYW3hCICIQzQdgQQgghhJDFAQVEFkHYT3QVpgLvpDdi7doNsm3bDrsSk2fyPaisdLzTLd0jHjXejYBADwnkUJQVTg0daqiYvKx2zoSlxv7VAqFgeyD8mgMxMjKsfw8M9F/1vgkhhBBCyMKAAiJDxOcEwPsAQz4mB8Iz6YEACFOCB6LIO/neVy945WhnkfzsmEN6B0YkELIvGUQGmtCVFU69hCVehxRMaJNxv70dDP6rBUIB3gd4HUpLy6SsrFzcbo+cP39W+vv7rnr/hBBCCCFk/sMk6gxii4bI/03pA1EQVYUJuN1GQEy+b2A8SmQU1IovZCdBwHmBkKeyIuwjNnm6yCvidYXFF3RqCBO2QznYq/0s3d2d4vF4pKamNiKQioqKpbf3inR1dUhlZdVVHYMQQgghhMx/6IHIEkZLBEKT4sATlQMB7I7UsR6IaFBpyezH67J09b++erLKkwHeCZNHgfeIw6UhTKaa01yAABkeHpLCwsIY7wr+LigolI6OyyzrSgghhBCyCKCAyBCJyppiFT8YJSAK4kKYgO2BSHxZRqOayMHDgNX/0pJiKXTHCoNCtxWTRxEIu7Qc7NUkUqMkLd7v9U7N2i4uLtH8iPb2y3PePyGEEEIIWRhQQGQNO4QpGGXre6OSqKMFRHESD8TQeJT3whWWkpISDSkqjgp5AgWukBR4JvcRCM+9mdzly21y/vyZSI5DImGklaM8Hrl8+ZIehxBCCCGE5C/MgcgytgfCMa0HorgA20wNN0JFJkOBy5Li4mJNvC72WNIbtR28E8VRXgwkXnstu3dD2dSIp2lpbT030Qnb0lKyyYAXYmhoQDo722XJkmWzOwghhBBCCFkw0AORIbBSH51EjfCjcNjS3g6GRAKitrZOKsoKE+7z/JWoBGy3JR6PVxOv49vHWVYoJgzKH8JYZNaVmLRsrA+9HwKaA4Fch2SYBHAIjtjkcUIIIYQQkk9QQGSBXxy8JN/7VVB+ccojwagIn/gyrqC5eb3s2Lpl2v05HSLLKoNqtCN0aGnF5E7fusWjYUTRidimF8RsKzENDw9r8zjkPoyNjemxpqOoqEhLvOJBCCGEEELyE4YwZYjoXIHv/NtJ/bdjyKWeA+By2OVVE1FQ4NUqS/AcJGJZlUOKvfaqP0KYVlaHxeF2SqHXKc31TunvD0tRAfZt95wYD4i4St0yNDSY0thDoaBcudKtVZsQvjQ+jtCnmWOf4BGB6MBxkGwNj0Syz0gIIYQQQhYmFBBZxj/hLHA57XyHREAUFHockW3jQYUllHuFgEApV3gytixBjoJ9OS0rHBPC1DEQlg11Hi3DinCkmTwJ586dkUuX2qSiolL/v6lpSUqfzS5D65BLly7qPqqra7TDNkUEIYQQQkj+wBCmDBOfD2BZtlfB5bRUBCQCRniyXhCgwAPvBVb33VJYWKSCAB2izfFwyFVNFerFABd6wuILF2gYUl9fz5T9IeQJ3aThecD7Ozs7ZHCwf6Kvgy0KElVfSkRpaak2lUPXavSGOHjwNeZEEEIIIYTkERQQGWImg9vthKch+ekv1hCkxHhdxgPh0dV9VECCZwEg5AjPVVeWS3PdpAtj0OfQJO7+/v4p++vr69VyrQhbgsiAcIA4gQiYrfcAYUxlZeVasQm9I/r7ezWsiRBCCCGE5AcMYcoAY76g7D89IsWO5D0REL3kdCY3zkuLcGlsr0I8KNOK9yKECcBgh/EPAgE7uRrPlUcVTRrxWVJZ4pSxsZGEHghUaIIIOXnyqCZNI/xofHxcCgomG8eduByU186HpLLYodWkti5zyeo65xSxhPwHPEpKSqW394q0t7dJWdmm6U4ZIYQQQghZzB4IJN8+/vjjctNNN8mOHTvkQx/6kFy8eDHp9j/5yU9kw4YNUx5tbW2yEPnhf5ySf31tUH60PyijyGBOwJjfkTQHAlSWJu+5gNAkeAaMdwAdqaNLryKsCU3mSgomDftRn/2e0dGplZjgcYBYaG+/NJE8HbL7UWifCfsYCEP6xYmgNrO72BuWy/1h+dmRgDz1ml9CYSt5KFZRieZEsDITIYQQQkh+kBEB8bWvfU1+8IMfyOc//3n54Q9/qILigx/8oK5sJ+LkyZOyd+9e+c///M+YR1NTkyxEfv76pcjfF7sSh++gwlKyHAhQVZY80dnjDGveg1n5t70EE52ugwGprKzU8KaSgsn9j/jsqk3o64BtRkaGpa3tgr4G7wPKtcITgVyK6uraKcccTdLEunPQkkt9U5veRZd2hUBBmBQhhBBCCFn4pF1AQCQ88cQTct9998mtt94qGzdulEceeUQ6OjrkZz/7WcL3tLS0qMehrq4u5pHv1Xum80DUJGkmBzyucExXaLtkqkuFGkREaWm5iouKEk9MCBPyGpAoDRFx6tQJuXDhvL4H1ZnwvitXujQJOtG4eoaTi4T+qA7Z8dgJ2E7dN5OpCSGEEEIWPmnPgThx4oSMjIzIDTfcEHmuvLxcNm/eLPv27ZN3vvOdCT0Qt912W1rH4UaW8jzA5UqcTL25CSLAnXSc9dWTYUnxFLpDUlRUGHlvcTEqMbk19AgGe2lpsb5WUeoVpyMgYcuhnbDhtRgbC0kg4FORgH8DgXEZHR2W6upqTcZOlvzdM5zc+D92KSTbV7iTvtdUZjp8+IDmVqxYsXJa8ZRPmET56RLmyeKAc4FEw/lADJwLZCHOhbQLCHgaQHz4UX19feS1aAYGBqSzs1P279+vYU99fX1yzTXXyCc+8QlZvXr1nMbgdDqkqqpE5gOlSTwJ1611S01NWSQROp7mlQgjOhX1DAx4h+xdWyClRT6pqiqPfMbSUq8UFxeqN6Gw0Cs1NRX6GrYp8o5o+BIepaVFMjIyKJbll6Iij/T3WxIIjIrDYUlZWUlMwnQ8HYOJczmMB6K11yGbliXO2ygq8orH45Te3i7p6+uWmppyWbFihSwmysuLcj0EMk/gXCDRcD4QA+cCWUhzIe0CAmVAQXSIDYBxCrEQz6lTtpGM8JYvfvGLmsz79a9/Xd7znvfI008/LbW1U+PxZwLlSgcHpyYL54KhQft8RFNdYomEQzIwMJZ0JT7gD0lNSVh6RpyytNwv62rHpbyyRpZWOaS3F2FIYenrG4nq/eCUoaEhFSSjo0F9LRSCtyIsIz6XjAcsGR5BnoMlFy+2a7L16OiYXLhwSUZGxsTjgQBJnOgQDFnS1hNbEQrVl851T4Y1HTjnk1U1yc+D0+mR0tIK6e3tkXPnLkhZ2TQb5xFYRcCNYHBwTEIoXUUWLZwLJBrOB2LgXCDzZS7g2Kl6P9IuIAoLCyO5EOZvgLh7JNTGs3v3bnnppZekqqoqEgLzla98RfMnfvzjH8sf/uEfzmkcweD8+BImqlCEqCY0lIMxD69BIhwOl9yyLiRXRpxS6hiRogKn1JQ7tEkcHsgriP6MXm+h+P3d2ocB+RB4zel0S5F78vhDY0i+LlAjHtv4/UEZHh6ZGIPdJyIRl/tCWrYVrKlzyvVr3VJe5JSwZcn/+6JPfEGR9v6wXBkKSXXJ9BMP+Rrd3d0qWiAqMZaxsVFZsmRZys3qFiK4EcyXOUlyC+cCiYbzgRg4F8hCmgtpD7IyoUtdXV0xz+P/GxoaEr4H8ffRxiOExrJlyzS0acGTwCZ3O+1O0tMZzMhXKC10yerasCZNI6F8MgnZzmeIBiVX0ccBHg1UYNLjuN1S5JkcAMKYIOrgJYIHAoe3/00+DXDMtt7JSbym3qXiATgdDtm5alKDvn4+cd+KaFBiFoIBOREQmUeOHJSWluP6HCGEEEIImf+kXUCg6hKSZl955ZXIc4ODg3Ls2DHZs2fPlO2ffPJJue6662L6E6Bz8fnz52Xt2rWy8EnggZihiRzA6xANxkNhBIQtIhxTciewso9t4IEw1asgJGIFhKUCA+FldrWmMq3AlKza1dG2oHzrBZ8cvDDZEG9pVeyU2bzUJYUTWuZsV1hDpab/XE79bBAQPT3ofD2qgubw4YORalAoJ0sIIYQQQhaJgIBx+t73vle+9KUvyX/8x39oVaaPfexj0tjYKHfccYcahwhhQa4DuPnmm9VI/uQnP6n5EIcPH5aPfOQj6pW4++67ZaGTqHKpSz0Q0wsIeCfgZcC5MaVQLSus/48kcZRkjT/vMMwLCibDxiAyCj2T3gNUYgLapbq8Qr0BPh8SqqdWfIIh/2JLUKI9aLVlDinyxnpNPC6HrG+0PwsioA5dnNkLgTEODg5oKVmETaGCFMq8dnV1atO5ffteShraRQghhBBCcktG6kShB8Rv//Zvy2c+8xl597vfrcbyt771LTWI29vb5cYbb5RnnnkmEvL07W9/Wz0Q2Pa///f/LmVlZfKd73xn2qpACxlUX02lxwVEgS0gsGrvVGMbjd4gFOJDmOB5gCFeUlIaeQ4iozjKAzE8HqtmMIbGxqYp+wIDY1OVz7I474Nh8xKXOCd0xcHW0IxeCFxXNLJDczmU+C0vr1RRAa8EelPgNXauJoQQQgiZn6Q9idoYpijDikc8yG1A34dotmzZos3nFgvwQMwUwmREgfE4wAsBrwAM66ampTFCQffpcmkDuZqa2hgPRHQlsOkavsWDpOh4llcnFhCVJU7Z0OSS45dD6oXoGAjLqtrknw9iqLKyWsdnqlDBG2JCmJDLgb8RYjVfMaFWycrwEkIIIYTkK7R+MoyVKAfCYYnbnYoHokAN1egQJgCREF/+taqqWtat2yi1tXWR52DclhQ4xeMSCYRE+kdSExBtvSFtDhdPU2VyhxXEBQQEONkekjG/JWsbXBrilPizxZb5hRekr69HPy8+d39/vzQ2LpH5ysWLrXL5clskzKyhoUlWrFilr+EzwIsynwUQIYQQQshcoYDIMCjVOncPBC4PBIRdsckIiESr3vBA1NfHVrnCdqj2VFEkcmVYZHDc0p4O7iRGPTh+OSgvnJiax3DtSpd6QpLRWDEpLtAfAg/s53ev80rVDKVdAT5fZaVdyhfVma5c6ZRQaENKoV7ZwuSjgLa2C+olCQYDKhjwN7psQzSgqhTCsXbs2C1VVZW5HjYhhBBCSFqZ/72yFziJeivAA5GKYYwcBrzfrlyEPg1GQEzNWUiECREqL5wcQ/fQ1PEMjVkyOBbWnhWJxMMtG92ye/X0WrO4wCFNlVMFxpOv+CWQQEQlAnkQ8D4gnAk5MZ2d7do/ZK7AsIeXAL0m0sGxY4fk0KHXZWCgXz0MPt+4hpnV1tbr3+3tl3Xc7e2XdJuenisa5oSQLEIIIYSQfIEeiAwTSlBNKJUqTEYAwOuALs7wQoRCwYk+D6ldNk2iLi6V2uJ+OSt2QvqpzlBMKFLfSFj+4VW/VouCUIinusQhm5akdrzrmz3yz6/7NQ8imgs9YWmuT92TgM+Hz3r+/Fk5e/a05nusXbtBk+tnQ3d3lxw9ekjDu+AdmCsQbm+88br09/dJMOhXUYNk9vr6xohHAuf6woWzKlhQmhZC6NKlCzI42CeWFZRrr71uzscnhBBCCJlP0AORg07UbodjSg7DdB4EPGCoYjUbf8eXcE0G3rNq1WpZVhnU3hMguikc2Hc2GCk1m8j7UFmSenfohgqnvHOHd0qviKNtoagmeKkRDIZkaGhQjfb29jZd+R8fH4u8jv2hDCzKvyYD70XSOcKLUDb4woVz6jk4e/aUNtCbCb/fJ62tZ/X4CKnCfuBNQAlaNOSLbgSI0rjoBh4Oh7QsLsrkIo9jYGBAent7Zd++V/TzzPY8EEIIIYTMN+iByDCJ2hl43XZp1plAeIxpKFdQ4JW+vpCubM+m8g+M2UKvW/MgekfsUq4mMds0l0tESQHG6ZDtK2KPFf3eRCypcsqSKlSPsuTvX/bL0Lgll/vDKlyW16TuhaisrNQGc3V19XpMJFgfOnRAli5dLkuWLNNVfiQyI2xo27YdarDjPOE92B7dzCEuiotL9Lk33nhNw4wuX7bDixAuhSTtZJ4g7KOl5YS0tp5TcQDx1tDQqDkP2B/2Gw2EHcrTRpcehscEoqGmpkpzItAgcePGLbofQgghhJCFCgVEhknUEM1pBVLyQKBSEao1QUSUlJRJMHhxViFMAFWC8J4Sry0g4BAZweI7jOzYQkgx3LO3QAo8U4UCukdjER1jg5CAcZ1IUCBn47pmtzx31I7/P39ldgICBr5piof94xxAEMAgxzEhKOxmhONy7NhhaWpaIsuXr1KhAOMfogIiA0IC1wBeA3gwkF9hexbOSVtbq2zZsj1htSQcBzkYJSUlKhpQItfkn5SVpZaDgmPh+MXFqKbl1AaK2CcFBCGEEEIWMhQQGSaYIISppCg1Qxqr2QhXQiUl5AGYkKbZCQivrrKXFEDI2KLl9fNBLdNaVeKQQHDq+HaudCUUDyAUwn4sNcKxKo+wHvSlSCQiVtQ4tcEcTsG57pDcsNY9bQWo6aioqJwQAv2aZ4AHzg9CiiBqADpZQyRgXBAVGCPEAcaG7bA9xEVFRVVkOzSzSyQg4KlANSjkTkzncZkJ817tyVFermNFGFV8Hw9CCCGEkIUCcyBy4IEoLylKqdMyVrth/GP1HSEz+H/8PRuDFu+HiIjuSG16PPSNWDIcV+RoVa1T9jZ7kn4WeDRgDNuhVB5dYUdCcSIQAmUStkf9It98wSc/O+wXfwLRMhOTwsmhIUgIR4LhjwRpCAIY5ngOVZAwHlRAij5X8FqgeR1KxUJI4PzDkO/ouCwnThyN5FLgM0KcoJISPtvViIdEHgk7F+O8ejWOHz+inhFCCCGEkIUEPRA5SKIuLynQ0JiZgPEKDwSM5+LiYiktLdXKP7MFCb9FnlHIiWm3qytzyI3rk4fnwOiFgMDKve0dgbhxRp5PBMq/XuqbTFg+2x2W0sKgvGldamFA8eC4vb1XdAXfGPcw9E0jOuSWQFRgjMmMf9NzAh4GeFAgHiA4qqtrVTicPHlMhQjyL9IJjotxI3cDx4Z4gQdk61Y7jCqVsDZCCCGEkFxDAZGDJGoYr4jXTwUIDSRC47F1q50sPJeV7/ICeAmSJz3A0N+1anrvBoxeOyQqNNGzwdJyq8gXwDESvRceiLds8mh3aiRTgzNddjjTXFb3YYAnCv+Bd6KsrCLS4XqmMrnGgwIxgs8FjwT6RSA3Al6Cqiq7qV26gZgbGxuR/v7eSAjVkSNv6Ph37Nil4yKEEEIImc9QQGSYYAIFASM31dXm7dt3xQiPua7al3rDUugRGU/S0wxN4GYymGFoo3IRkpNXr25WYxer6FittysTFSd834Ymlz5+etAvF3vDMuIT6Rm2pK0vpCVeUSr2zm0ecU3T6TpV43w22J6dEs2vQCgRQpkghuDByGQHbCRkI3wK+RUQEfB+4NSfPt2iVZrS5YmAOEIYF7wzzLkghBBCSLqggMgw4QRdmLMdqgKvAcRBY4VTqyEloqZ0+jHZ/QssNa43b94Web65eb3mHly8eCGpgDCsrHWqgAAvtgSkc8A+Nyj12tIRSrlhXTox5ViRRwEvBHIlMikeAAx6IwYxFxAqhSTvtrYLmoMBb1N9fYN6JfD3bIEIgkA5efJ4pBnh7t03zFpg5Qr06Ojs7JDGxqaUu64TQgghJHtQQGTZA7G0IrUSrunEzk+wZHlNYgGBng+FSaouRRt1ECLxIVQQJvX1TXLpUpsawclyIcC6Bpe8fCYowZBExIMBTexKChyCUVQUO6S8KLvnyK6SNPvwsHRgekjA6EcCN841xASeR34ERNt0YHt0x4aQQ2gZumaPjNhN70zfjmPHDkltbb0sX74yoacJ26HxHnIxcA3hvRgcHFRRtWbNWm2IB2FVW1sX8z5Us0KPi6amZdNe+1TBOTh16qSWu0WY1/r1mxKGdUFo4RxBiOXquhFCCCGLFQqIDBOK80Bcu2QspSZymfBArKpxyH+qlIjFmUKsPwxUu6fB1OTvmppaDcdBDsF0YVYoDbuyxilnuhJ7QZ55w46vKnCLvOeGxH0o8hmcO4gBhFQZLwLCqjZv3qphT4k8IzD8jxw5pOVt8TdClvA+NOLD3wBCoqOjXa5cQQ8PS69XdOlaCAdUhDJlb+GxGB5GlTBroht4rwQCQU2WR57JypWrVUhAPCB/A2ID261cuSYSKjXb/BFU1kISPDwP2B88L0g2HxkZ0cIBy5atUNEAoYcQura2i9pZHHNy167rI94vfHaMF+cLgqa7u0s9L9hHJnJaCCGEkMUIBUQW+0BsagiK12VpX4dsYjeTc0mB29Jk6aOXglpW1dBQkdywwgo04vRhhKILdOKmcU7tDg0DFUbedB6WJZXJBYTBFxRNuF5dl9lQovmG7Ymww4ywqg6jGBWijh07MvFcuTidbm1EB68EjHbkoKCELYQHun9jbsX3rkDODQz+gYE+OX/+jJaR3b37ejW+kQR/5swp9SLgPchnwXNY9YfYgKiBJwLXFeOBkQ/vBpLm8Tr2CVEJzwmMf/wNobN27QYdoyljnGhO4P2YXxj/uXNn1PhH1TGU2zUVvrq7OzXRHMntEDSmChi29XoLtbv3vn0vaaf2qqpa3R5iyiT1Q2zAu4MCBA0NTdOef3xunEPsf6GEexFCCCG5gAIiw9iN12yQI2y6KmcTGJAwyGAY7VpdILtWu2XUb8mP9vk0nAiiIhlILjY9GLAyngwYqFh9RlnU6UJKllVPNpcD9+z1yq9OB6VtIjfC0Dmw+AREPIj/xzmHgY05g9V0GModHZc0tOns2VO6cm97Feom8lSSr/5jVR7XE8nwBw/uk+3bd2s3bogH5GHgeIkMfiNqADwDEDXwPvj9AfVmYH5hbuG9SAiH2MQxkGgPwx/j2b5950Qo1WENpYKBjzAlGPzwfMBLhs8QPXZsg2Ojrwf2CSAYAEruYttQqEi9Y3genhTTKBDvwa6qqmpkeHhQTp06oa/hFEEchMMhPZ4Brx09elgFEr4rO3bsTthgkBBCCCEUEFntA+GICIjshzDBAIXRZCj2OuTeGwo0nClZ9SMYhDAOMWYYl9MJA6xYw2A8cuRgJO4+ERXFTnnrFo9cGQ7LjhVubTa3fbkroYCYr1wZCst/HAtoI77VdU41Sneuckt9uX1dgyFLRZLzKqtKGSPaGPYwbOFpgLHc0nJCV9fxmgkrmylEB+83XgEYyocPH1ABglLBJll5prmJY8B4RzhTtLAwxjYMd+wD4Www2jFGHA+hTngv+mzg2GgICOGAscPIn65nB7wOJr8iPoQOHovi4qm3seh8DIRd4ZgtLcdVbGB8eB9EDTws8DzAewZhhmuJc4yxX3PNzpiwMcxrfB6IKIwL5wBeF2yD0r8Q0QtBdBjPDwTX0NCQ5riY8seEEEJIKlBAZDEHArYZjMpsCwh4D7DqirCPaGYycGFImVCWZPkP0WDlF4YbVounM0jW1Lv0YVhS5RSvW8Qf1dC6c9ASX9CSAvf8i1t/fkI8gHPdttDpGPDLPdcVSM9QWJ47ikR5kd/aUyDlRVc/frsDt43LhaZ5XjXQca5RNWq2sf2Yf3gfPBsQARApswEGc7JKVWas8JzYuRi2YMGxkPuDVX+7Czs8BZVZqbKEz4u5iZAvfF6cL3hQDh58TcOyYPwb0WvnodjCADkYy5ev0AIBEEZIUocXCKV2IcYvXbog3d3degx4YJB3smpVs/4/xDauD74/86VBIMYI7xHCzRAWBjGFz4WQrw0bNif0MM7k1SKEELI4oYDIYhUm217PvoAAWBnFavNsgHEBg9CuUFQ5oxGBlWzEpWNleSZDMxp4QG7d6JGjl0Iy4rOkf9TSEKeDrUHZu2ZuDedA65WQNq3bumzSOzBXjl8OyuGLIRU+vRPiIRr01zh0ISgHL0x6eU5cDsne5vRfa1MN6WqMOhi28SFDmSxVC6+FIRdVkzCHo6tZQVDDk4PvBM4nzgPmuBk7zg+8EAiz6u3tFctCZ3Onzm0ktkNMQGiUlaFqle2lQz7IoUMHdB8QJtge+8Z5QOK5w2HJ8DAS0kWuXOnR7VEGGWIu00AInDt3Ws6fPxtJlkfIIcQdxOiBA/tUREx2dy9Wr+XJk0fV04JQzKVLl8nSpSvSUm2LEELIwoYCIoshTKh2hB/obFdhAjAWzGriTGBlEpVsYBghNn716rVqKM0EDEUYtjBU7OTrcCRGfiaMV+JiT0h+OlGN6UBrSB/N9U5583qPhl2lyoWekPzrIXs/V4YC6h2YK0fbgvJii+0e6T0X5SaJI1o8gP3ngtLeH5ZrVomsrpG0kg7DfzGvKuOzTzenS0tLNVwJHgdsB8+FCSOD6MCcjl+xh7fFfMeQawLDG98hrPAjNwQiJBSCd8qtyep2rohPQ6Vm6wWaDRgTSt4iUd0WNNURT5EJa4N35ujRQ1qBy4R44T4FLxIEA87XiRPHNKcFSfoYL5pHLlmyVMVh9KKInX/imLJ4gLyYVBcVCCGEzG8oIDIM4uHjcyCyXYUJmKoy0+UnGJB0ilVRxEhDeMBgSBWU20RSLjwtCIHCPiBCUjVWGyqmnhtUbeoe9Mvde7yRfhWvngnIpb6w3LDOow3y4oE3wwCPwbFLQfUA4Wogr33LUnvVeTrgDYEXw4iHeJD/cH2zW/79SECuDCcWZxjjpb4xFUDbltFwWkhgFd5UcwLGSJ6uVLHZNlpc4DsHEQHDua6uTvr7B9XoxvM9PT2aH4IqUZmo/IRjoMP5+fOnVRQkytGwxVS5elXwvYWIgLix811qIp/bDuPyaUgXQO4IqoChAhvCuvD54NU5fvywhnjBG4kKWM3N6zRM7OjRN/ScwhsTncCOfRqRApGB96QLLGRg/7hm0d93nBcjCAkhhMweCohFUIXJhI1g1RRGQvxqp51UORqppY+VR/zQo6LSbFdGYYiggg1izLEKizjyRMdMBpKqEzE4bsm3X/TJndd4pLTAIa+32gLhqdf88uvXemRplUvFWveQJTWlDmnvi03C/uXJWBHgdopsnKbz9amOkPz8eCBSLSqReHjbNtuzsrfZrdsijCmZk+e/WgIa4nT3bq8UzcKTQnJLujw9yMEw+U/4bqFcLIBhi0pUaPSHMskwpKPDveKBV8SIGnxnkRwOrx8Mbwh3HMc0+4OBDI8BPIIwyqfrao79RX9HixJ8XU04I0IVDfCg2P06hjXHAl4XiAhsFw736edE3gUECcbudKLfR6+sW7dBn8MY0aPEXqSwIrk9GAtEBp7H2FPNlcE9C8LmwoVW/RcNELGQgc+OcSOPC4ICHhfc93CusOhh+pcQQghJDQqIrIYw5SaJGuCHEz/GSA6NN+axOjg0NKAeEvxQezxuXSWE4JhLvDpWX/FA8y4YOe3tbZHk1VTYs8Yt+84GBY6at23zRBrMgX87FJCmytj9PH0gIG9aZ2m4EJKaCz0i/thooin84kRQuoYsefM695QqVIGgJS+2TBUPFUUOGRizNOH7hrWTX50VNS75/Rtdeq0HRi0pK3TIz08E5Gxcv4uhcUteOh2Q2zaz4g2xwQo4jH67gV6PGt7XXHNtpDIVjHKEH6H3Br6bqGLV2LhEjWwkQ0Ok4/0wxmEUI/Rq48atuhhw6dJFNe6RzzGdeLgasO/SUuR0dE0UWyiO8VpAICDUCSLAJPxjzGhcaHeu96oIse8/8H6gl4lbDX98VuyvuLhUVq9eo+V/E91DIJ6wSAFvA/JTzIKIHYrlUs8OFjQgciAe3G6XHhfbnTrVo4JmzZp1eq4rKioiVbaSgfdBnOB6QLDgM86lmAEhhCxkKCCyGMIE8YDV/VwICDtMoUJDDuLBDyAMfHhLRkb6Ncyhrq5BDZWrPeaqVWvUuJhNaMKOFS6pLXVIVYlDyoucsqQypI3lDO39U5f5f3Vq0sMAT0AqHLsU0uNsXhr7NTjZEYqpCAU2L3HJTRvccrkvLGVF9rjigRCpLrWNiJs3eMTlCOhzwz6Rtl5b0ZzvDqvQSFY6lyw+IBpgdMO4xSr9gQP71ViGiIdBDIMbq+jwXCK86OzZ0xNN9dA8r1pX/GEo19c36LYozwuBgX1hZT1T4sGA7zW8BYnua3gu2mMBsCgBY90uhRtUL4zdowPiYdI7i20gBOxO5wOyfv0mDYc0ixEQDMhHaW09G+kPgvfYXpxJEWOPsXCK8Q/sksa9cujQ63qPQrU53P8QdoVcERwD59CIA4gdeHVwHy0sLNZx4ZpgUQZ9VlasWKXiraurU4+PxRnkk+F4TU1LI3lihBCy0KGAyLIHAgIiuixnNrEN+KmGK35Q8QOLH0P8EKLjcbpig2EsmAZzqQoIGNcrayePD8MdpVMRnjRbfmevV3MUUNkpES+fDuqxSgocEQPj0MVJ98Wd2zxSWeyQimI7AX5pdWrnBbkat29B/w2Ehnjkh/81rOIDnhF4SZZWOcUftLQvBiF2CJGdM4DVdHgXYHhi1Ryr29F5F/GlVaPzLWzjPKjGcCarbMUz20URc38xBRYSjRPbmER3eFfgtTh92qMLG/BSoIwuQqaCwZB6XnBvM31rkmF62kSPG2ID3hAIB4g4CBYIEwgvnEfkrUDo4P8hHIx3A+LA3NNGRiBmBnSxBOOG18X2qCDPwq3vN14aXC90RUdSvdMZFofDG2muGI0RWTgWhA56lWB/WNyZL6WBCSGLFwqILOdAmAonucD+sbZX+uyVOHu1DgYHjPz+fnTnrdWEznSBH0CsuCFWO5UE7kRUlTg1d+A7/+mTsQnvgtsl8q6dXvnRPn/S9yGUqLrEIb+5yysdA2EJhOxEdggKAwz67/6XTzY0ubRkLP42wMhflYZu2PjM6Pb9kz57rL84gdhvO1/its0eWdfIRE4Sm4sA7wKMUIQyxYv5VAzkbPTXyCYQUfDCBAJBLUUL4C0oKSmLlOGdK3ivETLYJ4QchAO8DzDcEUKG/8c1wXYQc/HXBNcM20KImFAo423G/vEaxAnuv/B42FW54EVyqgDC50CYKbwUdpU+S8PP8FuB/iIQlJgPEBAQUNgG3hj7b3ozCSHZhwIiq52oHTnLgQB2Uyu7IzVW2PCDhh8gGLL4QcQPGPIW0l3nHcdoa7uozcRwfPy4zlZE4dwtr3FKS4ctyNY3uqSuzCkfurVAfvK6XxvPGZCfgHyIbcttw6LAIxGPBn60Gysc0jEQ65E42R6SwdHYnIVrlqfPsIcYQbgUqjUFo/Iz0NEa4VnbV7ikkt4IMgG+j/aqOjHAgMfDFHvIFHrPKCiYdTUo3NeTvcdOnrc9DPBkQGSgGl9JCRp8jsvg4KB6MZBQj8Rzu+CF3TUdyfG4J9u5GaJeEDv5/PKEyIBQQQhboz6wzdjYuOayTQobhJCmv8oXIWTxQgGR5RyIXLqe8aMEwwQ/XmZFDT8upitwY2OTlmRMNwgN2L59p4ZjYDUP4QiJut7OxA1rPVJTGhKPSyKr9gh3+o2dXvnpQZRSDcud27za2Xr7iuTGwW9c69UwpZfPxCY6tEeJClRyWlGTvmuF496+xSNPvjLVY3L8ckgf/8/1XooIQhYBuA/j9wDgN8EUq8D9GPcKhLHhvpno9wILQQDbwJsBwYBFICTiI98E+0DIFLaDlwKeFNx7ly9fJaOjw/o8Fouw/1w0dSSE5AcUEBnGlGuMzoHIFViBwkoWwgDMKhuEhFm9iu7Um26wbzzg4seP3lxA+dPtK6ZOWYgIlHIFqbjz8cO9aYlLG70FY50OEW/BHVvtuvTpBKFY6xqccqozwUFF5OfHAvLr13rF7br64+47G9AmfBBbN23wyNoGhkkRMt8xgiGVste4n0d7FZAPgvwMiBMsFqHyFJoH2jkUITlx4oguHiFXBLc2hK0uW7ZSQ6TssCyaA4RkivDE4gAw/+K5hZzPxDtGhglOKePqzHnJSLi+kciNHwz8yODHJtOVWgzwPGAc+LHD8eGNgKhAIuPVxjHPhgKPQ3tKXOwNy5nOkFZKMnkT/+0aT1qM+ERcv9YjvmBAyoscWtnpH1/1a3M7gDCsb77gk/pyh45jWbVTNjbNLr4bje8gjEzCuS8o8tzRgCaC15bNbu4hyRuldLNZMapnOCz/esgvLodD/tt2JLAnHrPm8ISRC8P4b0IA8l6iPQrxIUvRBgzAvbel5fjE70Kl9sNoaGjMu/yZ2WL/JiLc1773IEcQi2zMNZkf4N6PoiywGxDOaOY5coQQtofIimTAMxdtgyWyx0yuJrZF1TWEe5ueNdgeYhs2U8HEcVF+267shmIJIR0bmvGiyAK+V9geVeJQzt7+DtqV6RCiiNBC5D/B/kGeF/KvHA5EhMyu/1auoIDIMCENEbKBIZZrtWkS8uxkS5cMDtqu8nR2f53++FX6I4cfL6yA4QuFLyO+qNmO0V1W7dLHtSvtvhP4fdixwp1RoxTVnt6+ffIG9543FWh1JjSiM3QNWvpAB26MaWNT8q/p8LglR9pswYDmeBBEiZrf/dM+v6xvdMqtGz2R0AnDwGhYTneGZcRvaQld5IO8dt5O1IDQWVLp1I7aeD88QKbZH973X6fsDt8oW1s8UclqrkD8QOwg2R1hdU8f8Msd27zSUB77nRkPWNpAEJ8dTfyuWe5O+AMA79+wz+7LwR9/stiJ/+0xYgPGUX9/n3qGL1+u1epd8EzDmEGoE746aG4IwwgLPyhvG51EjvedOnVCvSZ2CJVdjhfGU7KFIZQYRm4HXkOfD2yL7y1CrUyYLfZnurPDmIOBlsr3eLbFOnAs/B6h6hnCwCAgMCZUAMNxkbuHc4KQMfxOlpdXSlMTerGwn082wPVE0QH0woFIgEGPeYJ5i/lhqqXBhsD/r169VoXF6OiYBAIw0gO6Dd5z+fIlvYa2J86epx5PgRZ6wRy/cqVTLl1q0/1g+8niB2YsYRXYeF9VVbW+bhdEgJ1n97HB3xCc2AavmSpqRpSi+lp3t12tzd6mW7q7OyIVOm0P4vVSVDT/wwspIDJMVBGmeVH7H70gMEExofFFwBcJPyTZLPeIpGr86OBGDdGA8aB0IsSF3Y8ilNVStyi5ijCfXADjFhWgYDyf7Z4a2vSL40HpG7Fk50q3ek2iCVuWPPOGX3pHUitviwT0lg6f5oAgTwSc7rSN9ujeGNEMjlkyOGY/B1GBsCiIDISTvXY+GOm5cf6KT8O+GiqcKmTixzr9uELy6tmADI/HPj/iE/n/9vvVe9Jcj/06ZNRvyeG2UEz/D/Ts2LLUpeLp5dMBGfWLVt8aC1gy5hdNvr9ji0c8SbqcE7KYQVgrjCEYy1hNhXFjG0GoIOWMhDvBELPL3ZZGmhNCCKD3CIxre1t75RaJ3XZvkiZdWYXwsA28In0vKmnBIMT2poAHKhbiN8E2siz9LcDiEv6GEQcjD4YhjEUIE1Qow9jxuuks3tnZruVyUeoW4zCVsSAG4PGGYYjfFxh4GDPGhPK4MBbtcdvCCMfr6enSv/E7idw9nAeII4fjgly8eF6bD6Icb64XBecruB6YHwirw/XH3+Z6YOUddgDOOQx6rMJ3dl7WVXuINSx0orAAzjOut73Kb9/3YRtEzz/YEpgHuHZY6T906IBeE3gF7HHY1Rfxr12cYDxyne1qlAMqZu0S2j59DXlD2CfmNL4bxj7C9hgz9t3V1an7hdg0+zNll5OVY54ObId9j4wMTuSmzn8clikqnkfgRtTbO5Kz47//oecjf5cXu2Vw1J7Iv36tW5ZWuWXv3jflbGyYpL/61S81JhY3PzSl2rx5mzZ8yxa4Gezf/5KOZenSFSpg4BXBjwDAFxxjy4dVY6z2FxWhJKc/Jh8mHqyqH2wNqgGOnhWvxCV4R3PrJreG9mAVPpkoeccOj/7bM2zJvxz0T2mMd89er3oS/uFV35TXrhb3RAfxpkqnih94MSAoTEEB4+HpGgzL4YvBKTkhtWUO9S6k2hAwFRAKduum5CIRXcKPXwpqmWAkz69OQ/neuc4FsjiYr/MhegUfhgwexrC3ywvDIMRNw+4aDqPM5GxEewpgjMEwN6uyZp+mlDmEBX4D7KIedsUpkyBujE8YmjDiTL6e/Zq9wgtvB5r34f/xm4L3wDiFZx2hJNHn2T6/dvNFe592aV27hG+Bjj9VDwfeh98q/AtRg0U5eC+wj5oaNDB0qXELr0r0/uzKW4nvK243DNUS6esbkWCixLw5gPGZPiQYL1bDMV5EHKQTGNMQCabLPOYA/kYneRwX4sGYmbYAk4gwM6IN1wDX26zw41zCe+Dz+fV1XGszN1IZj6ksmWqugRGpmMu5FoT9/T1yyy03S0lJVdrmwmyori7RBYBUoAcii2VcCzxeWbdufU7Hgy/Vhg2b9F98aUyJwGyC1aCiItwQcD426M0WN+/Dhw9O9KQo0R8CU5t9MQAvCPIjDCtrnfIPCSo2Ga9EMhBuBPFgvF3Ip/iDmwrkxZZgjHfhn1/3q1E/nXiAANF4U9/kvqM7gicD97yfvjFp/Rd6RH5tq0deOBFUYVBa6FDPRjwY8tZldj8OhGRhjMlAlSyIo1Q50R7SBygtEPX67Fljn+9A0JJ/fs0X+ZyoiPWO7fBcuPL+3oTfV2ceCHWSPqKNXhhTeEQbnR5P8tCK6N+SRKVwYUQbIxHgvo9HohDa6e7/2I8JgwLwLNiljxNXrgK4l0GomM8Sn9CaKtg/xAGECLwXWIm2f08RAnZRjVd8JnjUkVOC/0cfJDQaxLGXLl2m74eRPNdFMggmiBaEU+EzYRUeggy/nViQQygWFuUguuA9wvFRJhjHhwfJrKrjOuDzoP8TvEDxQDjiXOP3GAZ2MGjfk+38AL96DBA9gNexPzxnjHWcE/y2m07udmPFyXu2EZv41w6Ns/ur4P+xiIjQIgjE2Z6jRMUAZhIFppEnmR0UEBkmFFXGdcmSpRmtdJQqqLgE8KXHlxY3j2yCVSCIGCR0mx8Ok5jkcBToTQ0358UkIOKpLnHKW7d4YsKLkrGq1qklYq8MWSoY4kPlcHPcs9otHf3hSLgTkqt9QfvvkgKRd+6wqz9dGbK7ZJs8BzDmt7SaExYlYNRH99BAg7/zV0JysDWUMPcCwJPw9IHJzxEvHrBviAaIB/Nj0VTpkN+/sUA9A0isxj46B8LS1huWLctccn2zW5sAnu8O65hRwheCAl4PhDNhLKUFDukeCqtwiQZCAeFYCMPyuBzq+THiwfDCyYD87l5nJOwJP3zIEwFrG9BFXOSl00EZGAtLICh6vtAs0ISGzWfRgJwbhMshdAzeIRQNKC+a3+Mm+YG9An/1wtyIhdlgVx4sTOvvWPzvuTGQTbVBxLfbVa9sTwfCqxBnj/fa+YC2AAiHg1JXVy0jI/DE2KIK4D1YyYfhjtwTCAf0E4FogYAypdiNhwfeGggTvGavIjs0HAe/s/X1Dfo88l0wLju0xzbqYTzD628EBjwKGCu2tfdvh+5MJuLbzWjt4xVOGN+OiFjAa4mM9mgxYMQm/o0WnraoZA+chQBDmDIcwgQDJTAhIj72W+tl27r091nIBzANEVqFGwxuYrhpzgexlQxTax03PrPikmil5GrDFJCjcLYrpB4KNMRDDkI0iPW/a5cdjpQK8AD80z5fTHjQO7Z7Ul5tH/VZ8nprUJvioUs3xmXoHAyrcYqwJOQhpAKa/m1e6tLvSabm1bOHA3L+yuxdwUgaL/I4pLUnrGFlhroyCJPE17KhwiElXofUlDrV04KmhTDSTRNJMxd8/rDgthB9/jIJ7kGvngnG5I9Ee4hu3zx1DkBsDIxa2vUd8ytbY10szNcQJpJZ8HuBlXoTsgPDG3PB43GJ3w8j3Q4hszuV24tpdt6imSPIDynSfBATX29Co8zqPYRSqnmEOA48DUg6RgUg02wW48Iingntsv+1PVIkczCEiSQMYXLxi5cU3KDWr0dolVNdrZcuXZx1NY1sApew7Y6fXEkxbmH8QKBahN0x9+rGj/4N0T0cEHrzT6/6NUwIYU53bptdvwoYtTAWf3bErnaEbtuzCdVBpaUb1ycOeUO1JDwQ6vRvh/y67ZvWuuW5YwFNZobhjaICxgvyrp1e9RhkEpwb5D78V0tAPTSr65wqJhIlnl+3xi1r6p1aWhfn13Q9jyeZeACd6p2xYhLiEW6FnAp4Q/whv4z6wmqYw/NyxzaPLK/ObKgUvkfPHwvIuQRJ+gBiEiFnhZ6AigVN7G90yeutoYi3CFMMvVO2LUPHdId6eDDr4it6EUKmBwb45Kr9VDGJBVCDnQCePIwnXiTMZfUe+0aoER64V0y3IEZINBQQGSZ6YYk166cHXVMBQpjslZerz4OwKzFMrtBc3b7s1R5TIhTxmbhhI871woXzGiOKHwWIC7uMm3NCRKQPJE8jbKhv1NLQpbnc5CEY7n2TUwbGLKkvS/+cRCL47980+eP47uudGjKF3AN8HSAmir2zjz2eK1g5v33L5Dzaswalb0OaEwJNj9K6zfUu9ThgTNevdct/tsycWY5FGoQtQeC9ciYQCXGKB6FVPcNT9wcB99yRgPzWHngpMiOkYBCgL0i0ePC6RUVgU4VTXmwJyIUe+zX1SgVsL1V7f+x4IZRxvvBAyBuuYZHX/vwQG0ZIYMEEnxeeMdz7jGcM1bOQNL+82jkvqtERMl+Jvi9mO4w32ptByExQQGQRJiumht2XolDjPWGQI35zLglndn3z3okksbqrHhfiRrE6gwQ4CAP8i6oby5ev0trUeN0uX+jWz4C66YhXTXeUYHWpU6pLr96ozlZICoxIGK0AR4QBmkswj7Ytd+sjEVuX2b0u0JsDl+6m9ajNLVJT4pQ3LgY1TGvnKuRsTL7/rVu8csNaNLez5LVzITnZYYcKYTUfeRzJgLD6l4MBuWun96r7aMSDeffiyaAcuzwZtoSStrds8KgnCiD/4dWzwSmhcdMlrKO8rvkX+SVnOsPytms8Mu4XeeaQXytvGSBy4TX75UnbC2WOubLWNlJwvkZ8qNTlFF8A4RN22Gc6QCNE9C/BNYhvoojXUC2MHhRCCJlHAgJG1le+8hX5x3/8RxkaGpI9e/bIZz/7WVm+fHnC7fv6+uR//+//Lb/8pR0D/453vEM++clPqvGVT6QaV7bYgcFfXl6uBjhWYDCfkJAWXY8ZIG4T3opksZ5IHkPi2cBAX0w41GxDo3B81EaHmIFnBGNBDfBrrrk2sk1tbb0+j/mORLIlS5bJ8eNHVAS5XA6pqJi/+RxkKusbXdJQjvwllJV1xgiFZPMHngxIpLdsdsqN692agA0DGd4MdDuHIV1b6ZGAP6R9LZ563a+hTAgT+v5LPl2dRyI5mhsmAoY2EtoxHngJDrcFNSRsTZ1LDl0M6uu7VrtlVa1dYQ2lgKPFA/JN0AgwGnwOJLAjZ+NCT0jDvBDeBbsa48V7IKZ6h8NyvD2kSewYc7Rnta0vLN96IS4LfQKEi8Xnn/zroYBUlQRVMCTLlUF43upapyalp9oIEKIAQudUZ0jOdYciggXA44XrA3GDXeG6Oiau8w3r3MzvIISQ+SAgvva1r8kPfvADeeihh6SxsVH++q//Wj74wQ/K008/ndAld99992nM+Le//W1tHnL//fdrFYC//Mu/lHyCORCpg+oTiP3csuUaXSV8/fX9E92rkeDjVqMeTWOQNIYQokQGhikJi2Q1JJahogQ8AxAWECMzdRLFMVDHHKYGKlOgeZHpktrUtDRm26VLl6tHAuXosA1qgyMRHB4QvBdCAlWnmIC2cKgoTnytUjFmTfUmhOvcstEjN21Ah3NnTNLsr+/wyv/3mk9X8iEEYGgjnOiuXQ6pL7ebaSHpGZ3BkS/ROYjKJ3YoGMKNTH4dwrEMzx4KyLUr7f1cifIa3LbZo8Zyss8Tn2uTyOv15nX2+TAeNRzj348GJKrkfspEeykS0XolrA+AqmBvXu/WymQJu9QOW1pR7N+PBJL2DoFQQRhVzHtF1FOECl7oRg8vGY41U4iVEZAo/4vCAXVlzlk1TiSEkHwg7QIC1QWeeOIJ+bM/+zO59dZb9blHHnlEbrrpJvnZz34m73znO2O2P3DggLz66qvyzDPPSHNzsz73F3/xFyo4Pv7xj0tDQ4PkC1iJJqmxenWzhh2ZSkzoTNnRcVl6e336HCpTIIwIP+Yw1NGGHh4CiAqYBqhpjYoRTU3LVAi0tV2MtLNHB1OIiIqK6QXE4GC/hlJhnxAI27fvUmELIYCGPIm6uUazalWztLbis6yQ06fPSXd3Z6T+N1lcJApfRBjR27d75d8OBSJhTljZ//F+v6ypc0rXUHhKd24QX3I2GuwFyc/RwIuQTDzMBSOgEIZ0+2bRKlcGlLhFfgUEDzwAr59HiVw7X2dFjUteOBHQjuGzAWFIP9rn19AnGPgmBwmC6j9bAjP2A4kXXPHg/aZcMrTDugaX5sEgxAkfNTp3DRXRcExcJzSLRSlhw65VLtmx0q37u3AF+VKiHh3kGqFYwLJqu8AAygcTQshCJ+0C4sSJE1pF54Ybbog8h3CUzZs3y759+6YIiP3790tdXV1EPIC9e/fqj8Rrr70mb3/72+c0DnR2nG8UFrjn5bjmI6Wl6DI62dimuXltpP416mYjfGn58hVq3OOHuq/vioYrjY2N6Co/wpoKCwuktrZaqqsrZXR0WL0aBQVeqa2tlTNn+lR4oDZ2okZ6ppvm+vUb9f/hfcC1Ky8v1UcqLFu2VBob66W+vkrKyqrkyJHD6gVJd2I1WRiYePvouPu6cpfcc51TE53fuDCZuBxdxSkZ6NiNPaHE7Poml4ZCXZxIiDbctMGTNNcjHTQ3uOUur0ONZIiU6NV7j1vk1s2xIv0d13rlxOWQjnlwNCx15U6pKXPKxSuhSFng8mKHhnxBfJgqUPDQIFcE+TNNlS5p7w9FcjHiQYnht0wct7rU7otiC46wnOoIaQjYxiUu7SHy3BG74pYBwgBeCZPDAhBCVVwg4gvY5zoZ6C2CRyL6R0PaoNB0RUehASSVFxeGZHm1Q5oqpg/TshPbQ7KkyjURKkfy/d5AFiuOBRPynvZflo4OuzNkU1NTzPP19fWR16Lp7Oycsi3CnLDi3N7ePqcx4EuItvDzjarKknk5roUAzltz8wrNkzl37pyKzj17dqpH4OjRo2rc296BMhUEPT09smLFCmlsrNH3v/nNN2iVJCRWV1VVSXd3+8RKZkC7YsfT2zssDQ11smXL+rRUpVi5cqk4HCF5+eWXpaDAzVCmRUwBGnpEgVSvt24vkJu3WPKrk+Ny8LxfDWbDjlVeuWalVw3HoTFL+kZCsqzGLaWFkyFFZlX+WFtAznYGNITqzRsKpTxJGFY6aZ6lQ23X2qnPbYpLj9tVZm837rfkX14fldZuW1xBNKA3SjwQBSvr3HLD+sKk4URrSkTWRP3UbFkp0lDllZdP+TQfI2zZnoN44B0aSuAJgpcFX2MIi9kQ3RVdJCSvnbX/2rO2QFbXuTVv5lxXUB+r6t0qtM50GnEZkKXVLr3+21d6pYwNAPP63kAWHx7cWHThvWjxCQjkMoD4XAdU0kEr90TbJ8qLwPYw+OYCfkgHB0dlvjE8PC5ed+4a3OUDXm+xXt/Kylrx++FpKpb167fKoUMHZHwcDXe8sn795ki4UF+ffb49nhJ9AIQ0YTt4KewqTR4VCcaox+ujo+PS3LxRBgcTWA6zAKsIuBEMDmKel4jb7ZXe3n5NtLbzJewmPST/wcIGDASfL5C0cdh1a1yyqqZAfnnSLx6nXVK2sRI/KCGRsEh5gf0QKygTt9oYmuvwMEZI4m0WGndsdcuLJy3tpB7tAUDXdVRvQkUsk+geDgZkbOYKvBFKtImeO6ZR4sunUZI3lDDkCR6Q3as9sq7RpT0zEJqGxHZ0Q4cXyfaYWJoXgt4fq+ttj8Hzx/zSO0Oo1b7TPn1E03duapb5pd6QPl455dMmhQjrqipxyMYlk8ngwRAaAPK+kk/3BrI4CKDCg4ZQj8X0BMkWsFdy1kjONEdBLkR0oxSIgURVlbANto0H219NqEcuOvjNiOWYn+NaQFRX18m6dZtl6dIVkXPp8RTIjh275dVXX5LGxiapqqrVB0h2viEy7NKrV6S/v19/9JFzgbmIcqzIc8D/p+t64UYA0dLQsETOnDklDodLj4/vRHl5RVqOQRYGMBCmMxLQcO+3dk/Wul3sBgXs4Fs32qLoYm9Im/WhQhZyCoz4Ttc5QlduNB5E0jaEAfIZkMwN7wTyF9ZFh2khD8KypGAi+XppVfKcqnv22tcTno4DrUHpGrQ0B6aixCVvnPfH5FLMBgiWwTH7zb86FYzkfCBPBmFRCJVCzklVsUOaG1z6edBMEee0oji16lZk/twbyGLAitgM891eTLuAMOFIXV1dGkJiwP9v2LBhyvao0vTcc8/FPAcjDkYdwp7yCSZRXz1IVI5PVgZYyb/22t2a55AKSIKGQG1tPatJzajwhOpKKB1bX98ojY1LNIE63axZs1bzINBpG/kXyOHAsTGGRLkYhJBJ0LV7eRYqIqMXhacIHg7bk5AuEF51/VpPTPfhbcuc0j0Ir4JdOQsdy/FTsXW5W6tuwYOBpHosCrb3o7dFSEUBREgiTJJ9+4Al7QO2uIBwON1lGyP7z9mvL6tyaoWwsiL+LhFC5oGA2Lhxo5atfOWVVyICAqVZjx07Ju9973unbI8eEV/60pektbVVVq5cqc+hKhPYtWuX5BPswJpZZlvdCFWZICJQ2hWJ2TDk4bpbt27jlCpL6QJ9K1B2FtWgVq1aoyF8LS0ntMwsQq5mKi1LCMkvEHa0tMqlD0OyXiNLqhzaG2P3akvL/iK5GqIEPTySdUJPBvp3oP8IPF7blrll3UQndkLIJAgHhCMgEBQpKZxaUc+K+q4GQnao5eC4JV6X/V1FSkNLR0i9hdgMiwOwBbEt/r+xwqmPhVgKOu0CAvkMEAoQBdXV1bJ06VLtAwFPwx133KGrvL29vVJWVqbhS9u3b5edO3fKxz72MXnggQc0ERZN5+666668KuHqEHRZZcLbfAJf+muv3au5EC0tx+TixQtSX98kNTV2+FOmWLJkuYqXuroGLUeLBnSo+oR8DJSYTcc8wU0NIVL4F99J5FwsJGbb7I+QfGKmuY/Xoz0jqICF/h8wTtDRHD0v3rgQ0pAsNAlEBSgYMDBa0BjQVLcC3UOWPH88IK+etat2oa8FqmohTGymqkD4np7rDuvx4NFA8j4a+sFrguMkA6FfKPOLBn/z4XsOIxGfA00IYfChqWN9hVPOdYXFH7K0Bwk+K/qhoDM7FwOzCwoLwIuGXCN439BHBh45VEirLnFomWb0c0Hzy81LXHaflpDdpBNzEVcLjTBRhAKXDvlLeK7Q65CKIjscEpzpDEvvSFj33TFgqUDHd8GA7xN69CD3qLbMoZXc0A8HcwZV5xKV3U7h0+lYkMeEn/4CZ4Fs2eGTkgVQb8dhma5AaQQi4eGHH5Yf//jHMj4+HulEvWzZMmlra5Pbb79dvvjFL8rdd9+t26NizoMPPigvvviiJk/feeed8ulPf1r/ntvxw9Lbm7tk5fc/9PyU51wOS7728Zu4wjxPwZyFIV9ZWZW0s/VcQHUoVJBCMneyeEZ8BeGJ2L//Ze2ajVAmeFNwEzT/zgb0uMD3DvtBojZulRAQiYoVzEfQ+A8dvZEbEp1HFQ+S3XFu5oMBkgomZMU0kiOLm1zOB3QWR9I3emz4kiSdI7cDXy0YSzDCEGLVVOVUgw35E4Vuh5zpCsmpBJ4PGGDoJL5thVuKPOgSjl4cluZkoJQujm26kGN/yMuAkEHyOQQIckeXVjv1uDDu4aVBQjqsFTT8K0ggTvA+rPTiPSgP3D8S1tViGGUwEs3nMcYl9oXPj4R5dFhPdh7igaG3e7Vbw8rSde/J93sDGj3CeMe1u9wflpoSp16jaIGA3CacW9MnBXMUoq6tz857SvW0YD5h24FRayKbYGZMdf35knLwtj2Ncu/btuYkB6K6Gs16nbkTELlmPgoIt9OSr37sZsa5LzJSERAGhDJdvnxRm+J1draLZYU1pwNN8WYSP/BgQCxA/KAPC4xvJJSjbC26YF+50qX7gkBK9KNni5hRFe1Op0tDrNBrA3cHhFYhxMp05I5+PzwoyONA1/B0CC/koOCz2BW0ejXEDLko8WPG54LHBiFh8Z8pxqUcCKhHJx2leK+WfDcSyMKcDy3t6F8RjFlpnc9gBXlFjVPW1KNviEzkg4TVyEwFrBZPFLq5KjCGN61z675gHCMMBd6f+TgXIMJAJitz4b6LFXj0TEF+TuuVkAq0EZ8VEYvR4BbtcaIXjlPPH8J+YLfC+4PXIB7mAwjxg6diZNySPoiSaS5PWaHtmYCXAoL2bFdYv1dopLl5KQow2L1kIJSLJ4QS5i7OU++IvePSgrD8j7s3y5bmZfNeQGSuwxCZctNbIAulJIfdt5FYDeO5oaFJLlw4p0JiOk8Ebtp9fT1qIKPzNgz9FStWycaNWyOhUPBunDnTohWnENYEQRLtjcDrpus2Si3j5oHjQXAgPwRVqWCoQ/xCjOBfPB8KBTW/CY0isV+ICLw2l1U5iBCMHSJk40a7ylZLy3Edt9M5rGO2y++6NF8FwgDN/eyxD+jY8ZmQxwJRg/Niiwb0RwhrXxl8BkJILGhCuLbRKUfaQtrAD4nc8ADAOEdoyEw4Jhr2wQiEvQNDumc49RX92QL7GvkfeMyFROIBK9CIV68sRrCx3WUchh3CsFAtCyvkeA6dxWHoAoSuXOiZtIztkBjbOwIjGAZjVYlTDWjkqMAb0lRplx6+WjAe26ZwaKd3DBrjNOD8wzsEg3bUJ3LgQlDHd81yl4a+YaUfhvBcBAXurfAoBcIi57tD6jmCR6d/xIoYwantx+7kDuN58nPZ1zbR9YFRjvAxhNitqHWqRwuheTjfMNrx3lfOBNVAh/2rnrIih3o8cK6W19gllSGmIPRwfOyjYyCsxQkAvG2YA+gcj/mPuYxrGB16h/N5rjukYgKvN9fbggdnMv63b+8aWywgtM+wtCr2s6Gym6nQhnGODPfKqsYFEL9EAZF+kjl0jOuUkGTA+MZqOkCiNQxyeAMuXbqghjHC37BSVVhYrMYyVuptr0GhGtzDw4NqJG/atC1mrqFU7Nat2+XChfNy7NhhFQvo2g1j3e55MazCBZWs4HWAF2LJkmXq2UAp24sXW7Ws8qpVq1VgwLDv6rKbQq5cuUYFy4ED+zR0CuPBfsvKpg8/ivdg4D2rV6/VClsYC8a/cuVqFQjYJyqzQRBge2y7adNWqaurl9bW83L5cpt6X0ZGhlT0YDx+v+2NgScGAgh9QZBfMt13EJ/RlJSO97QQks/AuL1muVsf0cBgutxnJ2s3VjrVgD7Zjh4ZltQiH6AM8eMuKZ0wzPDzB2MJhhmMOYgKNOKDgWlWoZH38JZNHg1RQm8MxLMjhwJGYFWxU/eNeHd0IodxitVbrIWgrC6MeBiIycQJ9o3wKhhtEDNel52kPqK9Oqb+Nq9rsEvzIn8jlS7Q+H2HgfviycCUVXXsHXkdGCc6jx9uC6nhG7+IjNCn5TVOWVljG45YqXY4HeL2WnreMGYD9gUDFwY/9g3hAOFy7FJIbQqsYhvPEcQCtsO5NgbxTN3S8ZFhMO9Y6Za1Dc5IgvDAKEK7wtIzEtZGicgLwPERgoRrPBsPDgxxGO4w6pGTgDHj30SL69HeIYwNn6e53iXblru0Mlrsfh1yw9pYMYa5gnHimOZ6ppJPh8+LOQVvw3Tb4vyUFopsi/ueJAP7Qm5EKizEJGqGMKUZ3HA/+Jc/n/J8kceSxz9667wIpSDzM4QpETB+jx07oqv9EBjGyEUSNgSDLR6Wy8aNWyI3ymQ3QLzv3LnTE70v+rQjN3INSkpKZOfO61RoQJRopQhX4rsejnHw4H4VNDDWt2zZrgURIHICgaBcvHhehQ7Cj7BfhFMBeABAX1+fChS8F4IDY4VoQkUqlM9NBESOCUXCuDEGiAwzHpwbO6TJ7ueBErzR5wDekddee0XDsiCwbG8NGuZURrYzOTAYM84BtpkpdGyhhqyQ+cFimw/4rkIYwGC9mlAa7Aer1ghdMknj8B7AlJyut8XgGPp52MYsDHusWM91HFgthkDoHAhPNBRE0i3CdeSq0Q7nDlt4ZLOPGI6JFXh8nkRiazbUljpk63KXiiST0xBNAHkuEJsu0dAzNOCGxwBH7ei3ZCxg911J9N58p7+/R2655WYpKaliCNNigx4Ikk6QD7FlyzW6wo5EaKzWt7Vd0FwJhBkh7GnFitUpJRPDa7Fhw2YVIRAB7e2XNBQJTfVMk8eZ8hhwjM2br1HjPzr3AB4QAK8AxAc8EhA/8AagshpECwwmvAeGPLweBQVF6kkwHodkQDiYggrYf/x4IArwgKhKBELCcI5scePR7ygEB/ItIBJwbuHlgPdjx45dculSm5w6dULFWbTgV9e936fiB8eF6MBnARAd2FcqXpdcg8+BcDF4d/AZ7PLBDO8imQXfGXgZ0rEfeD3wmA3xoUMwlOcKVouRSB3/vYJAgffjRHtIE7rhucGKO7wtEC3wHMzUNHC2+RnwusDQhifH6NAiL1bT7SR0vL62waX7vdgTlo7BsMbzw3iHd8h4UvDe6E7viYAHAg8IG3idEOqDcC0cH8fCyj+Y6bcIXipPVGligyPu/8n8hgIizSTz55h4RUJmA+YMYv0NEBFY8YcBjhCnNWvWzXpewRhvbl6vQgKr9bPteYH3J6uQZsrFwiOCHA4jWI4cOagG/vr1m+TcuTOaIA6jHsZrNoDQgnhBaBIEDnJLEIoFrwM+C7wma9eu1zGuXLlKQ57g4YGXwjQZxANixvbSOPR5vA6RASHW2Xl5ItkcpXMLxe126bGSleXNdqlaO1wNFbrGVHziukNI2cnqJRERSQiZPfgu2zX+RbavcOsjnj1r3JpgDHGBJGPE6MOYxm3AF3SIPxBWIx4Vo8Ym8k8QxuOYMPBhvBd5HLJrtVvfi9Av/IvFGYSQIbwKIUwQDPEhPwB5BPGgMeHhiyE9rhEQ6Ga+tNqlHcyRWwERZMLDaMcQA0OY0kwwFJY//OtfTHm+vFDk0Y/elpMxkYUbwpQMrJ7DCLyaykco9QqvRDb6kyDMCPkc+PExt5xc/xBdudKtBj+MfyMCDBAYR48e0gRt9AWB1wFCCEZ2V1enio6mpqWybNkKFXIAwqij47IUFhZoaepwOKShKfic2LfDgb8tnQejo2P6Ot472waIcwHnHOFf8DRA5MGLBDGFMZw5c1rOnz+jcwmJ+DBG7HkFkTSgnwHPQQwhpyTX1y1fWGwhTCS1uQD7Bd8xJELbCerZ6x+FHAyTX0ByQz9DmBYv04UwEZIuoo3duZLNcJvonIr5YoBO53lBFSwY0whlWr9+oyZgI4EbnorTp0+qIQ0BEf1ZmpvXyZo1a/Xv1tazGrqFsKje3iuaiI5zUF1dIQMDQ7JkyVIVcKhqZYdguTNbXnF4SMXD5s3bdEwIzwIY07p1G3QunDx5TJxOhHeF1FMBQwbjxPY4F62t5/Sz2OFyTvXIIGck1V4l80U4EjKfMd8PNADMNhQOZDZQQGQxhIkQsnBAbgaMZ5MfYIx8hGfN9OO/alVzjBiBkQ3PxLZtm+TChXYpLi7T1f833nhdPRoIr8pEgQUY7fAiwPuDcSNkLdGYly9fqeFbGD6S7eGtgNcFlbDM58f5wPMIeULoE8Zv/t8WGXZeCJoh4j6oIR0ul3pZED6F5+HVqKioYEldQghZ4FBAZElAcNGNkIVHOpKLESKGnBOEs8GghliAaxrP2wnyr2svD+RTwOCf7pjReRMwyhGCBW+CyW/A3yY/xfY8DKu3AMeBkEkG9hmdoI4yvvEgTwIP8xr2jwpXhw8f1IR5fB6E4iDkCw0A4WFB2BcS6JEPUlFRpaFsKAVsiwrklLg018T2Yjj1MyCsDucpPol9tkniOCcYA/JtjJfHhGkhDwSv47g4X/SKEELI7KCASDPJmqfTA0EIiQcGLCpRHTnyhq7Qw6hF0nx0yJdpnIfVeyTPw4iHsW4EgzHeEVaFsrymDC+SvBHLum7dpmnFw1wxJXg3bNik+SIYC5LiEcpl8jowRuSFwLsBDwbCohAq1dbWquFg+GzYhz1uh34GfFZ8BvQrwd94r51DMv1NFOcOHhGIFuR34Jxg39gnclngWcEDr+McQ8Cg+hbEBMZNEUEIIalDAZFmGMJECJkNMHCvvXaPVnGCsd3b26s5CUh2xgo/yszanb9DUltbr8YwVtJROhc9QNBID8Y7emm0tMA4v6iGN/YLYx6J3pmkrq5BduzYrWOM955A7KBSmAEOBQgmeClg4MOARy4JwqDgOUDjQggBJJV2d3dpyWIY+Mi9QPPDeI8E3oNSushPMSVp8UDeCs4fRIHJ98A5wfnE+UOlLdyrcYzjxw/LlStdum94RHDu8Z5kggLHwZjg1cD4jdcInhOEfZlxmX8xNuzLvEYIIfkABUSaoYAghMwWVHdC3gRyFC5caNXGfBATCPtBYjZyEWC4wgiFEYywG9PoDsa4AXkOEBkIXVq2bHlWKjzBOEYORKpA/MSHSEFExAPhgwc8KeigjlAoGPd2roZdRhehSjgP9fVNKgBMKd5Exj+MfCTORyfPNzUtUbHS3n5ZyyNDUECQoEcGzh08JHbBAlSicqlgwAPHwGsQK/g8ECVomojzbntWEEJlRXqYQARCJOH6GaGlfQPQWpkQQhYgFBBph1WYCCFzA0YrVs+rq6vl/PlzsmLFykiHbrO4P52xjlXxRInSCxnkSyCH48SJo5rT0N8PbwW6m1ep8EDuxtWEaOG95v3IVYG34/TpFhUF8ETAOwGRgP4fpaXl2r8E3h4IAIgX/Gt7SXrk0qWLE8KiQYUDtm9oaNTXUJnLeEyQCwKR4/G4JBgMSWFhiW4fXVLZdIWHcIH3BOFg+H0pK6vISMI9IYTMBgqINJOsmjcFBCEkVeBFwIPYYNV+9+7r1ZCGkY6Ve3hYMhEWhO7ve/bUqIcDRj68EuPjo5pXgpCxRCWUES6FB6pZmYaD0V4Qk3wOTwTyPdB3o6DAKw0NtdLTMyAXL16QsbERTfAGyM8wQCzBY4HEdIwBAgeiAqKGfTkIIbmCAiLNsAoTIYRkBhjSWP3PNCahG5hGgamSzDsAUQHBU1ODR22kyWRv77B6mRAC1dXVruIA3g88h5Atn88vdXV1Oh6Ik8uXL2llKfzb2dkRyTvBv0bc2InoyLvwqtcD78tkrxFCyOKDd5RsNZLTnpKEEELIJDD04fXAAx4MgCaEyIuJ9y5AK5hmhRAY7e2XVJTgd+fy5YsR7wTEBCpO4e3w1sCzgefg2TB/I1cjHWWKCSGLEwqINMMQJkIIIVcDwpVmAp6RaO9IY2OThlshlwK5NFeudEt3d6dWtkKOBv5GKBXCqCA47BwPe1t6Jwghs4V3jTTDKkyEEEKyDSpDrVixKvL/pifHZFPBIRUQ8ELAK3HmzCnp6urQXA8sfSHhmx4JslCIbqqZSfBdsb8z9vHivyOWZenDbqQZ1rLNfn9AwuGgeDwQ6IX6HEIUzRIzcpmwLUIR8TeEPDyJ0UUUFgIUEOkmWQgTBQQhhJAcAMPHeCtMjgb6cUBgoPneuXNn1COB36ni4lLNmcB2xjOB95vw3HijDYaQ3bQQPTSKoipMeac18EyVqeimibkEIV8o0QsPDcZu8kdQ+QsGnt2fBHksXq12RrIHrgOuC3J7cJ3sOeaMlFe2jfjwRO5Spc5FO/cnqHPZLq3s10pmk0Y6BIFT/x/zFsIA+8ccwHP42w7/c0z0iHHqviAmLAv9Zez5gf0CM9UxN7C/oqJK/V6gAhv2vWTJ0ki4IUo6Q8xDXNgloIcnGola4vUunLk1P765eUQ4WRJ1tgdCCCGEJME02UMJW+RfoKP5qVMntNs5DCsYTzCYYPCYMrbGmML/w/CHgYZeHBAn0BcwjJC7gffAeIIhZRvfds8LY4gZsYLO5NgHtjMrvQDbY9vJ94YioVYmITyROLE7jY9P7A/7cui/MOrM2KPBZ4LBCGPR7XapeMI2xkjFmGpqUJFrSI1AjBd/o3cIyukutBXjXGD6nWAVHtcDhj6uIaqWGS8CHnYzyLBecyNYbUEXnLg+br0GKN9cVVWjHjYUHejv79cCA3gfjHD0i7HnBq63U+cpvG8w3lHiGiICIhDY+/br67je2BZd7PE8GldiTjQ2LlVDH/MC48H70SemowMFD5yyevVa3R+OiTFi7pnqaOh6j7wkvI7vWKI5i8+JOYv+Pgg1RPhifX29jIwEZL5DAZElWGqPEELIfASGGYwrGGcItTAhGTBoYAShFwZWcVHlCUYbDDoYXDCM8D4YUfiN6+i4pIYS3osE7wsXzuvqKgQJXh8dxSqyXxPA8R4YZDCcYPhh/2jyh+1gwGPlFoYVwGotjDwYn1hxhoFnh314VWCgBK7tIYAIKNFjGrGBzwZRBIGD4xnPCPYFYYHPhH4iMEjRYwWvYV/o54HPDmMQ48C2eB65JK2t56Svr0dKSsp0DNgOxupMXpdcYgzkTGLCcmC8w9C2RYNt/8CAn1zRd2u+jseDUB6ICIkk9+Pa4Vzi2hUVlehzmCu1tbVSWFgcMdYBesCYuQowL9vb21SIlJSU6DExDyBmIXLjRYsRN8gXwhzH9YegNnMjHrzfsGbNOh33dN4ozNFETTKjwTggOkzTT1Rnw2ekgFiEJKvCxBwIQggh85loQyaa1aubI38jz8J4C+K7fqObumHt2g26emvCSGDkwZBDSBA6q8O4MiQKj4LRjpVePAcjH6FWtrHokN7eXhUoMNzxOjqRY+UWq8c4JpoNer2FKgJg9EE4YBzoOI7jwvOCfWF7rGbHeyZMWIs5F9HnBL1AkJh+5sxJ7RaPfcPgg7EMz4rpAQIhAwPWFkQ+FUUwhs1zV4sJqbLDe4IxXhm8BvFnwnfsJoT2+cU1MyvwZvvZjscWCvZx7eMEJ16xDXPMDQgvXGdcOxjZ6PSO64lzDQGB/BuMBeIP4rCurkHPGYRZV1enXreVK9dMVAxLbKQbIWDAdmj0GE9075bo7Y1XCs0e8Zhk5kaNHuYLUUBkrQ8EFQQhhJAFjjE+U8u7iG26l6yLeiIDFkY+Hgb0zjDAy4F+IDDMYbTGl7yFgW9jHx8hKCBdzRmxv82br1FjFcIAhi/GAuMX40AoF0JpEO5kQlswRgie4eFhNT5N3L0dKmZJQYFJrI09H9g/nsd5x98weO3Ed5lY7Ycx69UcAQBjHWIGxrl5Hiv1OCcQPBATWOWHEDMJwHboDgx92yTEc9gHxB7EAFb0TQKwvXLv0PdjLPDcQCzY3h80O3TpeGGQx6/OR1/PZP1c4ClIJALI/IMCIs1YSQq5OuepW5MQQghZaMDYzWX52ejEdACBAGMawADGqjrCsBBHj1AciCcY/mgAiJAthG5h5d7uKu6S8XGE/MBYR68O216A4Q6RAaN8eNjOKUFXdOwLnhB4ZeBRgTCCMIAHBCFW8BLhgX2bHBPsE31GIEIGB/v1NYgD/D/CfuDVsUWCLSyMyIFwsfNRCiOJvwg3wpjhTWE39MULBUS6SVrGlV8wQgghJN8x4iK+izmMbTxgoCPmHjH7MMy9XiTfumRkxC9DQ8OR7ZHY3tCwROP5IT7gJYAgwH7hXVi2bEVkWyMqNmzYHJfcPekxghCA0MEjmqamJdLTY+ecQFxA8MDjE+2VICQezow0YzsUpwJ3HyGEEEIWNxAYiPk3IHG2qgrx/yNSWVkTed50JgcQHqkwl8pQCDWKzQEgZGZYgyzdsA8EIYQQQgjJYyggshPBRA8EIYQQQgjJCyggslWFiS4IQgghhBCSB1BAZK0PBAUEIYQQQghZ+FBAZMkDwRAmQgghhBCSD1BAZCsHgvqBEEIIIYTkARQQ6SZZCBMVBCGEEEIIyQMoILIVwsQcCEIIIYQQkgdQQGQphIlVmAghhBBCSD5AAZFmWIWJEEIIIYTkMxQQaYZVmAghhBBCSD5DAZElKCAIIYQQQkg+QAGRZqwkWRAMYSKEEEIIIfmAO9079Pl88tBDD8m//du/yfj4uNx2221y//33S3V1ddL3fP3rX5dHH310yvMnT56UhQZDmAghhBBCSD6TdgHxwAMPyP79++XLX/6yeL1e+dznPif33XeffO9730v6HgiFd73rXfKJT3xCFjrJBASrMBFCCCGEkHwgrQKis7NTnnrqKfnGN74hu3fv1ucefvhhufPOO+XAgQNy7bXXJnxfS0uL3HPPPVJXVyf5WoXJRQ8EIYQQQgjJA9KaA/Haa6/pv9dff33kudWrV0tDQ4Ps27cv4Xv8fr+cP39e1qxZI/ncB4I5EIQQQgghJB9IuweiqqpKCgoKYp6vr6+Xjo6OhO85ffq0hEIhefbZZ+ULX/iC5lDs2bNHw5nwvrniducmPzyZp8HpcuRsTCR3uFzOmH/J4oVzgUTD+UAMnAtkIc6FWQmItrY2uf3225O+/j//5//UvId4ICggDJKFL4GioiJ57LHHpKenR8Oe3ve+92k4VGFhocwWJCxXVZVILijpG0/4fHGRN2djIrmnvLwo10Mg8wTOBRIN5wMxcC6QhTQXZiUgEIr0zDPPJH39hRde0JCkeCAeIBAScdddd8nNN98cU6Vp3bp1+tzzzz8vb3/722W2hMOWDA6OSi4YGhpL+LzfH5S+vpGsj4fkFqwi4EYwODgmoVA418MhOYRzgUTD+UAMnAtkvswFHDtV78esBITH45Hm5uZpqyn19/eriIj2RHR1dan4SEZ8iVeELlVWViYNe0qFYDA3X8Kkx7VyNyaSe3Aj4PUngHOBRMP5QAycC2QhzYW0Blnt2rVLwuFwJJkanDt3TnMjkNeQiEceeUTe9ra3xVQvQqhUX1+frF27VvKnjGu2R0IIIYQQQsg8FxDwMrzjHe+Qz3zmM/LKK6/IoUOH5OMf/7js3btXduzYodvAO9Hd3R0Jdfq1X/s1uXTpkvaPgNhAtaaPfOQjsnPnTrnpppskb6owsYwrIYQQQgjJA9Ke5v35z39ebrjhBvmTP/kT+cAHPqDlWR9//PHI6+gHceONN+q/YOvWrfJ3f/d3Gv5099136/s2bdqkvSQWYvO1ZH0gWMaVEEIIIYTkAw4rmcW7wGPHentzk7B89Hyv/M0PD055/ndvXiJve9PGnIyJ5A6U7kX1LSTQz/d4RpJZOBdINJwPxMC5QObLXKiuLkk5iXr+F5pdaCTNgeCpJoQQQgghCx9atdkKYWIOBCGEEEIIyQMoILKWRJ3lgRBCCCGEEJIBaNammWQZJUyiJoQQQggh+QAFRJZCmBZiRSlCCCGEEELioYDIVggTBQQhhBBCCMkDKCDSTbIQJp5pQgghhBCSB9CsTTNsJEcIIYQQQvIZCogshTAxB4IQQgghhOQDFBBpZkVDqSTSCuwDQQghhBBC8gEKiDRTW1EkD/3RDVOep4AghBBCCCH5AAVEBmiqKRG3K1YwMAeCEEIIIYTkAxQQGSI+54FVmAghhBBCSD5AszZLAsLBU00IIYQQQvIAWrUZwhV3Zl3MgSCEEEIIIXkABUS2PBAUEIQQQgghJA+ggMhWDgSTqAkhhBBCSB5AAZEh4gUDk6gJIYQQQkg+QLM2Q8QLBqeDp5oQQgghhCx8aNVmCJZxJYQQQggh+QjN2myFMDEHghBCCCGE5AEUEBkivugSqzARQgghhJB8gAIiQzjjBIOLMUyEEEIIISQPoFWbrT4QDGEihBBCCCF5AAVEhojPeWAnakIIIYQQkg9QQGSI+IgleiAIIYQQQkg+QAGRIeIFAz0QhBBCCCEkH6CAyBDMgSCEEEIIIfkIBUSGiPc40ANBCCGEEELyAQqIDBHvcIgv60oIIYQQQshChAIiazkQPNWEEEIIIWThQ6s2Q8R7HOiBIIQQQggh+QAFRIaI1wtOeiAIIYQQQkgeQKs2SyFMdEAQQgghhJB8gAIiS52o6YEghBBCCCH5AK3aDBGf88A2EIQQQgghJB+ggMgQsYLBwjO5GwwhhBBCCCFpggIiSyFMhBBCCCGE5AMUEFkQEI4ESdWEEEIIIYQsRCggMkS8XqCAIIQQQggh+UBGBcRnP/tZ+fM///MZt2tra5MPf/jDsnPnTrnxxhvl0UcflVAolMmhZRw2jiOEEEIIIflIRgREOByWhx9+WJ588skZtw0EAvKBD3xA//7hD38oDzzwgPz93/+9fPWrX5WFTIzHweGgB4IQQgghhOQF7nTv8MyZM3L//fdLa2urLFmyZMbtn332Wbl8+bL8wz/8g1RUVMj69eulp6dH/uqv/kr+6I/+SLxer+RDDgQhhBBCCCH5QNoFxMsvvyzNzc3qQfjoRz864/b79++XLVu2qHgwXH/99TI8PCzHjx+X7du3z2kcbnfu0jtcLqfE943L5XiI5HQuRP9LFi+cCyQazgdi4FwgC3EupF1A3HvvvbPavqOjQxobG2Oeq6+v13/b29vnJCCQf1BVVSK5JDpkCX/mejwkt5SXF+V6CGSewLlAouF8IAbOBbKQ5sKsBASSnW+//fakr7/00ktSXV09qwGMj49LeXl5zHMFBQX6r8/nk7kQDlsyODgqOfVAxOU89PWN5Gw8RHI6F3AjGBwck1AonOvhkBzCuUCi4XwgBs4FMl/mAo6dqvdjVgKioaFBnnnmmaSvR4chpUphYaH4/f6Y54xwKC4ulrkSDIbnVRWmXI+H5BbcCDgHCOBcINFwPhAD5wJZSHNhVgLC4/FofkM6QfhSS0tLzHNdXV0RwbJQYdElQgghhBCSj+Q8S2PPnj1y7NgxTZqOTsQuKSmRjRs3ykIlPoSJEEIIIYSQfCDrAgLhSt3d3ZGwpbe+9a1SV1enFZtOnDghzz33nPaQeP/7379gS7jGhzBRShBCCCGEkHwh6wLiwIED2m0a/5qE6W9+85vafO6ee+6RBx98UN7znvfIH//xH8tChg4IQgghhBCSj6S9jGs03/3ud6c8d91118nJkydjnlu5cqU88cQTkk8whIkQQgghhOQjOc+BWAxQSxBCCCGEkHyBAiILUD8QQgghhJB8gQIiQ1hWrkdACCGEEEJI+qGAyBCWRCkIuiAIIYQQQkieQAGRBagfCCGEEEJIvkABkSEYwkQIIYQQQvIRCohsQBcEIYQQQgjJEyggsoCDCoIQQgghhOQJFBAZwmIMEyGEEEIIyUMoIDIE5QMhhBBCCMlHKCAyRXQVV0YwEUIIIYSQPIECIkMwgokQQgghhOQjFBBZaCRHBwQhhBBCCMkXKCCyARUEIYQQQgjJEyggMkVMCBMVBCGEEEIIyQ8oILKgHygfCCGEEEJIvkABkQVYhYkQQgghhOQLFBAZgo3kCCGEEEJIPkIBkSEoHwghhBBCSD5CAZEp2EiOEEIIIYTkIRQQWegDwTRqQgghhBCSL1BAZIjoFAjKB0IIIYQQki9QQBBCCCGEEEJShgIiG9AFQQghhBBC8gQKiCyUcaV+IIQQQggh+QIFBCGEEEIIISRlKCCykUTNOq6EEEIIISRPoIAghBBCCCGEpAwFRIbYtbE+8veWFSU5HQshhBBCCCHpggIiQ9xx/Sq5cWuNNNcG5ZZtlbkeDiGEEEIIIWmBAiJDuJwOuetNS2TvyqB43a5cD4cQQgghhJC0QAGRYZhATQghhBBC8gkKiAyLB+gHighCCCGEEJIvUEBkHIoHQgghhBCSP1BAZNwDQQFBCCGEEELyBwoIQgghhBBCSMpQQGQQ5j8QQgghhJB8gwIio9jigSKCEEIIIYTkCxQQGYQ5EIQQQgghJN9wZ3Lnn/3sZ8Xv98tDDz007XZf//rX5dFHH53y/MmTJ2XhQwFBCCGEEELyh4wIiHA4rILgySeflN/8zd+ccXsIhXe9613yiU98QvKxDwQhhBBCCCH5QtoFxJkzZ+T++++X1tZWWbJkSUrvaWlpkXvuuUfq6uok/6CCIIQQQggh+UPacyBefvllaW5uln/5l3+RZcuWzbg9QpzOnz8va9askXyDHghCCCGEEJJvpN0Dce+9985q+9OnT0soFJJnn31WvvCFL4jP55M9e/ZoOFN9ff2cx+F25y4/3OVyRv51Op36by7HQ2RezAWyuOFcINFwPhAD5wJZiHNhVgKira1Nbr/99qSvv/TSS1JdXT2rASB8CRQVFcljjz0mPT098vDDD8v73vc+eeqpp6SwsFBmi9PpkKqqEsk1paWFUljokbKywnkxHpI7ysuLcj0EMk/gXCDRcD4QA+cCWUhzYVYCoqGhQZ555pmkr1dUVMx6AHfddZfcfPPNMcJj3bp1+tzzzz8vb3/722e9z3DYksHBUckVUI64+MPD4+LzBWVoaFz6+kZyNh4iOZ8Lg4NjEgqFcz0ckkM4F0g0nA/EwLlA5stcwLFT9X7MSkB4PB7Nb0g38V4LhC5VVlZKR0fHnPcZDOb+SxgKWWJZdlWq+TAekjtwI+AcIIBzgUTD+UAMnAtkIc2FnAdZPfLII/K2t71NLFjaUaFSfX19snbtWlnYmM/ETGpCCCGEEJIfZF1AoOpSd3e3/gt+7dd+TS5duiQPPPCAnDt3Tvbt2ycf+chHZOfOnXLTTTfJQoedqAkhhBBCSD6RdQFx4MABufHGG/VfsHXrVvm7v/s7bSZ39913y5/8yZ/Ipk2b5Bvf+AaNb0IIIYQQQhZDJ2rDd7/73SnPXXfddSoWornhhhv0kY9QBBFCCCGEkHwi5zkQhBBCCCGEkIUDBUSGoQeCEEIIIYTkExQQGQb6gSKCEEIIIYTkCxQQGcQuTUvxQAghhBBC8gcKiAxD7wMhhBBCCMknKCAIIYQQQgghKUMBkWHogSCEEEIIIfkEBUQGqaiolMbGJdLQ0JjroRBCCCGEEDL/G8ktdpxOp6xYsSrXwyCEEEIIISRt0ANBCCGEEEIISRkKCEIIIYQQQkjKUEAQQgghhBBCUoYCghBCCCGEEJIyFBCEEEIIIYSQlKGAIIQQQgghhKQMBQQhhBBCCCEkZSggCCGEEEIIISlDAUEIIYQQQghJGQoIQgghhBBCSMpQQBBCCCGEEEJShgKCEEIIIYQQkjIOy7IsyTPwkcLh3H4sl8spoVA4p2Mg8wPOBWLgXCDRcD4QA+cCmQ9zwel0iMPhWLwCghBCCCGEEJIZGMJECCGEEEIISRkKCEIIIYQQQkjKUEAQQgghhBBCUoYCghBCCCGEEJIyFBCEEEIIIYSQlKGAIIQQQgghhKQMBQQhhBBCCCEkZSggCCGEEEIIISlDAUEIIYQQQghJGQoIQgghhBBCSMpQQBBCCCGEEEJShgKCEEIIIYQQkjIUEIQQQgghhJCUoYBII+FwWB5//HG56aabZMeOHfKhD31ILl68mOthkQzQ398vn/3sZ+Xmm2+WnTt3yrvf/W7Zv39/5PWXXnpJ7r77btm+fbvceeed8tOf/jTm/T6fTx588EG54YYb5Nprr5U//dM/ld7e3hx8EpJOzp07p9fzxz/+ceS548ePy3vf+169J9x2223yne98J+Y9vG/kF0899ZS8/e1vl23btsk73vEO+dd//dfIa21tbfLhD39Y7xk33nijPProoxIKhWLe//3vf19uv/12ueaaa+Q973mPHDt2LAefgqSDYDAojz32mLzlLW/R+8K9994rBw8ejLzOe8Pi4G//9m/l937v92KeS8e1n2kfGcciaePLX/6ydd1111k///nPrePHj1vvf//7rTvuuMPy+Xy5HhpJM3/wB39gvfOd77T27dtnnT171nrwwQeta665xjpz5ox1+vRpa9u2bdbDDz+sf3/zm9+0Nm/ebP3qV7+KvP/P//zPrbe+9a36/jfeeMO66667rHvvvTenn4lcHX6/37r77rut9evXWz/60Y/0ud7eXr0nfPrTn9a58E//9E86N/CvgfeN/OGpp57S7/r3vvc9q7W11fra175mbdy40Xr99dd1fuC6/uEf/qF18uRJ69///d+tvXv3Wo899ljk/T/+8Y/1PvLP//zP1qlTp6xPfOITuk1PT09OPxeZG48//rj15je/2XrxxRet8+fPW/fff7+1a9cuq7Ozk/eGRcL3vvc9vQe8973vjTyXjmufyj4yDQVEmsBFvfbaa63vf//7kecGBgb0x+Dpp5/O6dhIesEPAYzE/fv3R54Lh8MqCB599FHrf/2v/2X99m//dsx7Pv7xj+sNAHR0dOgN5Re/+EXkdYgQ7BOGBlmY/M3f/I31vve9L0ZAfOMb37BuvPFGKxAIxGyHHwLA+0b+gHvAW97yFuuhhx6KeR7fe8wDXM+tW7da/f39kdd++MMfWjt37owYBZgXf/VXfxV5HfPmlltu0feThcdv/MZvWF/84hcj/z80NKT3h2effZb3hjyno6PD+vCHP2zt2LHDuvPOO2MERDqu/Uz7yAYMYUoTJ06ckJGREQ1JMZSXl8vmzZtl3759OR0bSS9VVVXyf/7P/9EQBYPD4dDH4OCghjJFzwNw/fXXy2uvvQbBrv+a5wyrV6+WhoYGzpUFCq7bk08+KQ899FDM85gLe/fuFbfbHXkO1/38+fNy5coV3jfyLHzt0qVL8uu//usxz3/rW9/SsCXMhS1btkhFRUXMXBgeHtZQhJ6eHp0X0XMB82b37t2cCwuUmpoa+fnPf66hawhVwz3C6/XKxo0beW/Ic44ePSoej0d+8pOfaChzNOm49jPtIxtQQKSJjo4O/bepqSnm+fr6+shrJD/AF/mWW27RHwLDs88+K62trRqviOvd2Ng4ZR6MjY1JX1+fdHZ2qggpKCiYsg3nysIDovGTn/ykfOYzn5ny/U82F0B7ezvvG3kmIMDo6Kh84AMf0B//3/md35Hnn39en+dcWHzcf//9akQipwULTo888ojGta9YsYLzIc+57bbb5Mtf/rIsX758ymvpuPYz7SMbUECkCRiHINqoBDASkTBL8pfXX39dPv3pT8sdd9wht956q4yPj0+ZB+b//X6/zpX41wHnysLkgQce0ATJ+JVnkGguGOGIa837Rv4ATwL41Kc+Je985zvliSeekDe/+c3yx3/8x1pUgXNh8XH69GkpKyuTr371q+p9QGGNP/uzP1OPE+fD4mU8Ddd+pn1kg0nfB7kqCgsLIwai+dtcyKKiohyOjGSS5557Tn8QUFXlS1/6UuRLjHkQjfl/zAXMj/jXAefKwqy4A1fy008/nfD1RNfa3NyLi4t538gjsNIM4H34zd/8Tf1706ZNWkXp//7f/zuruRC/DefCwgOrwKiu9+1vf1vD0AC8EBAVWJnmvWHxUpiGaz/TPrIBPRBpwriaurq6Yp7H/yO2neQf3/ve9+QjH/mIluj7xje+EVH/mAuJ5gG+1FiNgtsRZWDjv/ycKwuPH/3oRxq7Ds8TvBB4gM997nPywQ9+UK91orkAcK1538gfzPVav359zPNr167VGHjOhcXFG2+8IYFAICZXDiAeHuGunA+Ll8Y0XPuZ9pENKCDSBJKiSktL5ZVXXomJjcbq0549e3I6NpJ+fvCDH8jnP/95rev98MMPx7gSsdr06quvxmz/8ssvq5fC6XTKrl27tMazSaY28dPIjeBcWVjA6/TMM8+oJ8I8wH333Sdf+MIX9HriOkfX+sdcQNI8Eix538gfkCBdUlKihmM0LS0tGvOO64nrakKdzFzAezAPMB8wL6LnAvoIwMPFubDwMPHpJ0+enDIfVq1axXvDImZPGq79TPvIClmr97QIQN1/1Ox+7rnnYur2ov43yR9QcnXLli3W//gf/8Pq6uqKeQwODlotLS36+l//9V9rfeZvfetbU/pAoKzrbbfdZr388suRPhDRZd7IwiW6jOuVK1esPXv2WJ/61Ke0rj+eR61u1Ps38L6RP3z1q1/V8osotRjdBwLf8/HxcS31/IEPfECvs+kDgXrvhieffFJLNWJ+mD4QqPXOPhALj1AoZL373e/WEp4vvfSSde7cOeuRRx6xNm3aZB08eJD3hkXEpz71qZjf93Rc+1T2kWkoINJIMBjUGt7XX3+91v790Ic+ZF28eDHXwyJp5utf/7oaiYke+DKDF154QRvNoe47fkB++tOfxuxjZGREmwrt3r1bHxAUaAxD8ktAAAjEe+65R+cC+gR897vfjdme94384oknntDFASwioA8AhEJ0Dxk0ocQPPWq4o28MDM1o0Hjy5ptvViHxnve8xzp27FgOPgVJB+j58cADD1i33nqrCsvf/d3ftV555ZXI67w3LE4Bka5rP9M+Mo0D/8mOr4MQQgghhBCy0GEOBCGEEEIIISRlKCAIIYQQQgghKUMBQQghhBBCCEkZCghCCCGEEEJIylBAEEIIIYQQQlKGAoIQQgghhBCSMhQQhBBCCCGEkJShgCCEEEIIIYSkDAUEIYQQQgghJGUoIAghhBBCCCEpQwFBCCGEEEIIkVT5/wFFaCL1Y4buSwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n = 1000\n", "\n", "uniforms = rng.uniform(size=(n))\n", "uniforms = -4 + 8 * uniforms\n", "\n", "m_n = np.cumsum(uniforms) / np.arange(1, n + 1)\n", "\n", "m = 2\n", "\n", "s_n = np.sqrt(np.cumsum((uniforms - m) ** 2) / np.arange(0, n))\n", "\n", "\n", "m_n_low = m_n - 1.96 * s_n / np.sqrt(n)\n", "m_n_high = m_n + 1.96 * s_n / np.sqrt(n)\n", "\n", "fig, ax = plt.subplots(figsize=(8, 4), layout='tight')\n", "ax.plot(m_n, label='m_n', lw=2)\n", "ax.fill_between(np.arange(n), m_n_low, m_n_high, color='gray', alpha=0.5, label='95% CI')\n" ] }, { "cell_type": "markdown", "id": "6c165f34-971d-43d3-92c9-5bc41ef1a5e6", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: LFGN loi de Cauchy\n", "\n", "Reprendre rapidement l'exemple précédent en remplaçant la loi uniforme par la loi de Cauchy. On obtient des réalisations de la loi de Cauchy en utilisant la méthode `standard_cauchy` de l'objet `rng`. Répliquer plusieurs fois le tracé (avec l'axe des ordonnées restreint à $[-10,10]$) pour différentes valeurs de $n=100\\,000$. Qu'en pensez-vous? " ] }, { "cell_type": "code", "execution_count": null, "id": "7289a19c-477a-4291-acda-ac8345bb7d42", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "68ecfdc4-7625-41e1-8bcb-9a2e178ba238", "metadata": {}, "source": [ "## Illustration du TCL \n", "\n", "On veut illustrer la répartition de l'erreur renormalisée $\\displaystyle \\varepsilon_n = \\sqrt{n} \\Bigl(\\frac{m_n - m}{\\sigma_n}\\Bigr)$ pour différentes valeurs de $n$. Lorsque $n$ est grand cette erreur renormalisée est proche de la loi normale cenrée réduite, c'est ce qu'on veut vérifier numériquement. \n", "Pour illustrer cette répartition, il est nécessaire de répliquer un grand nombre de fois l'erreur c'est à dire de considérer un échantillon $(\\varepsilon_n^{(j)})_{j=1,\\dots,M}$ de taille $M$ et de constuire l'histogramme de cet échantillon.\n", "\n", "**Attention:** en pratique il n'est pas nécessaire de répliquer $M$ fois l'estimateur $m_n$ pour approcher $m$. L'estimateur de la variance $v_n$ suffit pour donner la zone de confiance autour de $m_n$. C'est une information importante donnée par le TCL." ] }, { "cell_type": "markdown", "id": "31f6167f-b1a7-48a5-a766-21e39afbd9a0", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: TCL loi uniforme \n", "\n", "Dans le cas de la loi uniforme sur $[-4, 8]$ vérifier la répartition de l'erreur renormalisée $\\varepsilon_n$ pour $n = 10$ puis $n = 1\\,000$ à partir d'un échantillon de taille $M = 100\\,000$." ] }, { "cell_type": "code", "execution_count": null, "id": "092b36ef-6d3a-43ed-a1f0-e1cb3d4ddcd6", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "0a5b46ab-d911-4bd4-a94d-cea1cf217889", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "aa218f54-374f-4d2e-8773-8eff7162ab37", "metadata": {}, "source": [ "## Un premier exemple d'estimateur de Monte Carlo\n", "\n", "On va mettre en oeuvre un estimateur de Monte Carlo pour calculer\n", "\n", "$$\n", " I(\\beta) = \\mathbf{E}[\\exp(\\beta G)] \\quad \n", " \\text{où $G \\sim \\mathcal{N}(0,1)$ et $\\beta \\in \\mathbf{R}$}. \n", "$$\n", "\n", "La valeur exacte $I(\\beta) = \\exp(\\beta^2/2)$ est connue mais cet exemple permet d'illustrer l'importance des bornes de l'intervalle de confiance (et donc de l'estimation de la variance) dans une méthode de Monte Carlo. La seule valeur moyenne $I_n = \\frac{1}{n} \\sum_{k=1}^n X_k$ n'est pas suffisante pour déterminer $I$." ] }, { "cell_type": "markdown", "id": "227b7a6d-1c9c-41c8-a15d-235c3a309ee7", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: fonction `monte_carlo`\n", "\n", "Ecrire une fonction `monte_carlo(sample, proba=0.95)` qui à partir d'un échantillon `sample` de réalisation indépendantes $(X_k)_{k=1,\\dots,n}$ renvoie un tuple qui contient: \n", "\n", "- la moyenne de l'estimateur Monte Carlo de $I = \\mathbf{E}[X]$,\n", "- l'estimateur de la variance asymptotique apparaissant dans le TCL,\n", "- les bornes inférieures et supérieures de l'intervale de confiance de niveau de probabilité `proba`." ] }, { "cell_type": "code", "execution_count": 30, "id": "2f51f3ad-eb43-44a5-80ef-5bb5592115a3", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [ "def monte_carlo(sample, proba = 0.95):\n", " n = sample.shape[0]\n", " m_n = np.mean(sample, axis=0)\n", " s_n = np.std(sample, axis=0, ddof=1)\n", " alpha = 1 - proba\n", " quantile = stats.norm.ppf(1 - alpha / 2)\n", " ci_size = quantile * s_n / np.sqrt(n)\n", " stats.t.ppf((1 + proba) / 2, n - 1) * s_n / np.sqrt(n)\n", " return {'mean': m_n, 'std' : s_n , 'lower' : m_n - ci_size, 'upper': m_n + ci_size}" ] }, { "cell_type": "code", "execution_count": 33, "id": "fe42ea47", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'mean': np.float64(1.0247101435825938),\n", " 'std': np.float64(0.20850911066050495),\n", " 'lower': np.float64(1.0117868524848315),\n", " 'upper': np.float64(1.0376334346803562)}" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beta = 0.2\n", "sample = np.exp(beta * stats.norm.rvs(size=(1000)))\n", "monte_carlo(sample)" ] }, { "cell_type": "code", "execution_count": 62, "id": "4b72c09a", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0.20.51.02.03.05.0
mean1.0189451.1506381.5881495.22553964.10526019262.613472
std0.1961840.6015822.13876413.632066594.367456281052.583878
lower1.0067861.1133531.4555904.38063027.2666601843.113942
upper1.0311041.1879241.7207096.070447100.94386036682.113002
\n", "
" ], "text/plain": [ " 0.2 0.5 1.0 2.0 3.0 5.0\n", "mean 1.018945 1.150638 1.588149 5.225539 64.105260 19262.613472\n", "std 0.196184 0.601582 2.138764 13.632066 594.367456 281052.583878\n", "lower 1.006786 1.113353 1.455590 4.380630 27.266660 1843.113942\n", "upper 1.031104 1.187924 1.720709 6.070447 100.943860 36682.113002" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "result = {beta:monte_carlo(np.exp(beta * stats.norm.rvs(size=(1000)))) for beta in [0.2, 0.5, 1, 2, 3, 5]}\n", "\n", "import pandas as pd\n", "pd.DataFrame(result)" ] }, { "cell_type": "markdown", "id": "3f1c949d-63c4-48f2-9523-29445c96a71e", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: premier exemple\n", "\n", "En utilisant la fonction `monte_carlo`, reproduire le tableau suivant où chaque ligne représente un résultat pour une valeur de $\\beta \\in \\{0.2, 0.5, 1, 2, 3, 5\\}$: \n", "\n", "- la première colonne est la valeur moyenne $I_n$,\n", "- la deuxième colonne l'estimateur de la variance,\n", "- les colonnes 3 et 4 sont les bornes inférieures et supérieurs de l'IC à 95%,\n", "- la colonne 5 contient la valeur exacte $\\mathbf{E}[\\exp(\\beta G)] = \\exp(\\beta^2/2)$.\n", "\n", "Ce tableau est obtenu pour $n = 1\\,000\\,000$. Comment interpréter ce tableau? " ] }, { "cell_type": "code", "execution_count": 6, "id": "69aa7e19-9d81-4ce9-9ccf-a93221f56bed", "metadata": { "tags": [ "question" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meanvarlowhighexact
0.21.0205514.245536e-021.0201471.0209551.020201
0.51.1340273.650418e-011.1328431.1352121.133148
1.01.6510604.676691e+001.6468211.6552981.648721
2.07.4026852.379946e+037.3070687.4983017.389056
3.087.9150758.558333e+0682.18127393.64887790.017131
5.0121963.8256196.439313e+1472228.169429171699.481809268337.286521
\n", "
" ], "text/plain": [ " mean var low high exact\n", "0.2 1.020551 4.245536e-02 1.020147 1.020955 1.020201\n", "0.5 1.134027 3.650418e-01 1.132843 1.135212 1.133148\n", "1.0 1.651060 4.676691e+00 1.646821 1.655298 1.648721\n", "2.0 7.402685 2.379946e+03 7.307068 7.498301 7.389056\n", "3.0 87.915075 8.558333e+06 82.181273 93.648877 90.017131\n", "5.0 121963.825619 6.439313e+14 72228.169429 171699.481809 268337.286521" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "df = pd.read_pickle(\"data/first_df.pkl\")\n", "df" ] }, { "cell_type": "code", "execution_count": null, "id": "d0796346-e25d-4afb-9d58-fff5bf7cb7a1", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "da8bbc68-33fb-47fb-9816-de3b4f2d6cbc", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "048d17c1-5a98-4d72-aaae-5c7370165e87", "metadata": {}, "source": [ "## Option panier: un exemple multidimensionnel\n", "\n", "On considère $d \\ge 2$ actifs financiers dont la loi à l'instant $T > 0$ est modélisée par une loi log-normale c'est à dire \n", "\\begin{equation*}\n", " \\forall i \\in \\{1,\\dots,d\\}, \\quad\n", " S^i_T = S^i_0 \\exp\\Bigl( \\bigl(r-\\frac{\\sigma_i^2}{2}\\bigr) T + \\sigma_i \\sqrt{T} \\tilde G_i \\Bigr)\n", "\\end{equation*}\n", "où le vecteur $(\\tilde G_1,\\dots, \\tilde G_d)$ est gaussien centré de matrice de covariance $\\Sigma$ et les constantes $r > 0$, $\\sigma_i > 0$ sont fixées. Il s'agit d'actifs financiers $(S^i_t)_{t \\in [0,T]}$, $1 \\le i \\le d$, modélisés par un processus de Black-Scholes multidimensionnel. On introduit la matrice $L$ triangulaire inférieure obtenue par la décomposition de Cholesky de la matrice $\\Sigma = L L^\\top$. \n", "\n", "A l'aide de cette matrice $L$, on définit la fonction $\\Phi:\\mathbf{R}^d \\to \\mathbf{R}^d$ telle que \n", "\\begin{equation*}\n", " (S^1_T, \\dots, S^d_T) = \\Phi(G_1, \\dots, G_d) \\quad \\text{ou encore} \\quad S^i_T = \\Phi_i(G_1, \\dots, G_d)\n", "\\end{equation*}\n", "où $(G_1, \\dots, G_d) \\sim \\mathcal{N}(0, I_d)$ (l'égalité précédente est à considérer en loi).\n", "\n", "On s'intéresse au prix d'une option européenne (aussi appelé produit dérivé européen) sur le panier de ces $d$ actifs financiers, c'est à dire qu'on veut calculer \n", "\\begin{equation*}\n", " \\mathbf{E} \\bigl[ X \\bigr] %\\quad \\text{avec} \\quad g(x) = (x-K)_+ \\quad \\text{ou} \\quad g(x) = (K-x)_+ \n", " \\quad \\text{avec} \\quad \n", " X = \\biggl(\\frac{1}{d} \\sum_{i=1}^d S^i_T - K\\biggr)_+.\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "id": "cce85a35-c93f-4b1e-9369-12ee4a6272f1", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: initialisation \n", "\n", "Définir les paramètres globaux $d = 10$, $T = 1$, $r = 0.01$, $S^i_0 =100$ (pour tous les actifs), $\\sigma_i = i / (2d)$ (on dit que certains actifs sont plus volatiles que d'autres) et la matrice de corrélation $\\Sigma$ définie par $\\Sigma_{i,i} = 1$ et $\\Sigma_{i,j} = \\rho \\in [0,1]$ pour $i \\neq j$, avec $\\rho = 0.2$.\n", "\n", "Initialiser la matrice $L$ en utilisant la fonction `np.linalg.cholesky`." ] }, { "cell_type": "code", "execution_count": null, "id": "8e3bd2ae-8e53-468b-9adf-b816bf4d5f0b", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "8427fa3d-b0d1-402e-a8dc-b1e32ae620ec", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: simulation d'un échantillon d'actifs\n", "\n", "Définir la fonction python `phi` qui transforme le vecteur $(G_1, \\dots, G_d)$ en un vecteur $(S_T^1,\\dots, S_T^d)$ (tous les paramètres sont des variables globales pour simplifier l'écriture du code). L'appel suivant doit fonctionner \n", "```\n", "G = rng.standard_normal(size=d)\n", "phi(G)\n", "```\n", "Si on veut implémenter un estimateur Monte Carlo il faut travailler avec des échantillons _i.i.d._ $(S^{(j)}_T)_{j=1,\\dots,n}$ où $S^{(j)}_T = \\big(S_T^{(j),1}, \\dots, S_T^{(j),d}\\big) \\in \\mathbf{R}^d$. Modifier votre fonction `phi` pour création un tel échantillon à partir de l'appel suivant: \n", "```\n", "sample_G = rng.standard_normal(size=(d, n))\n", "phi(sample_G)\n", "```\n", "(il faut utiliser la technique du broadcasting en `numpy`, c'est très important à connaitre en pratique)." ] }, { "cell_type": "code", "execution_count": null, "id": "356d0d99-d09d-485f-b526-ddcdb2bff272", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "5d8ff622-af28-4b12-aba9-92a1b3b422fa", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: estimateur Monte Carlo \n", "\n", "Définir une fonction $\\psi: \\mathbf{R}^d \\times \\mathbf{R}_+ \\to \\mathbf{R}_+$ telle que\n", "\\begin{equation*}\n", " \\psi(G_1, \\dots, G_d, K) = \n", " \\biggl(\\frac{1}{d} \\sum_{i=1}^d \\Phi_i(G_1, \\dots, G_d) - K\\biggr)_+\n", "\\end{equation*}\n", "dans une fonction `python` appelée `psi`. Cette fonction doit fonctionner avec un échantillon $(G^{(j)}_1, \\dots, G^{(j)}_d)_{j=1,\\dots,n}$. \n", "Ecrire et programmer l'estimateur de Monte Carlo pour estimer la quantité $\\mathbf{E}[X] = \\mathbf{E}[\\psi(G_1, \\dots, G_d, K)]$ où $(G_1, \\dots, G_d) \\sim \\mathcal{N}(0, I_d)$. \n", "\n", "Pour différentes valeur de $K \\in \\{80,90,100,110,120\\}$ et $n = 100\\,000$ vous devez obtenir le tableau suivant: " ] }, { "cell_type": "code", "execution_count": 10, "id": "fbf3346e-91e5-4cbe-9c48-599ed101d177", "metadata": { "tags": [ "question" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meanvarlowerupper
8021.394471228.31877221.30081821.488123
9012.860460181.18794712.77703212.943889
1006.655165111.5537496.5897026.720627
1102.99865054.1329852.9530493.044252
1201.20415821.9762781.1751021.233213
\n", "
" ], "text/plain": [ " mean var lower upper\n", "80 21.394471 228.318772 21.300818 21.488123\n", "90 12.860460 181.187947 12.777032 12.943889\n", "100 6.655165 111.553749 6.589702 6.720627\n", "110 2.998650 54.132985 2.953049 3.044252\n", "120 1.204158 21.976278 1.175102 1.233213" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "df = pd.read_pickle(\"data/basket_mc.pkl\")\n", "df" ] }, { "cell_type": "code", "execution_count": null, "id": "c5919e9a-e6e2-47f0-b8cd-71f38413b8e7", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "469f3249-caf3-46a0-8fea-671a2887a908", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "question" ] }, "source": [ "### Question: monte-carlo adaptatif \n", "\n", "Ecrire une fonction `monte_carlo_adaptive` pour calculer le prix à une précision $\\epsilon > 0$ fixée (telle que la taille de l'IC à un niveau de confiance donné soit plus petite que $\\epsilon$)." ] }, { "cell_type": "code", "execution_count": null, "id": "ecc8a015-6811-4d51-aecb-65a043e4ee56", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "9416f9fa-c4e0-4597-8d9e-1dd3a872653f", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: variables antithétiques \n", "\n", "Sur le même modèle que précédemment, implémenter la méthode de Monte Carlo avec réduction de variance par variables antithétiques c'est à dire basée sur la représentation: \n", "\\begin{equation*}\n", " \\mathbf{E}[X] = \\mathbf{E} \\Big[ \\frac{1}{2} \\bigl( \\psi(G_1, \\dots, G_d, K) + \\psi(-G_1, \\dots, -G_d, K) \\bigr) \\Big]\n", "\\end{equation*}\n", "Calculer le ratio de variance (variance de la méthode naïve divisée par variance par variables antithétiques) pour les différentes valeurs de $K$. \n", "Que signifie ce ratio de variance? \n" ] }, { "cell_type": "code", "execution_count": null, "id": "02560fe6-d34d-4054-ac36-170b6c68cbfd", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "806f8959-f831-43da-ae57-cbb36e2ee66a", "metadata": {}, "source": [ "## Option panier: une variable de contrôle \n", "\n", "Dans le cas de la dimension 1 ($d=1$), le prix est donnée par une formule fermée, on appelle cette formule la formule de Black-Scholes. Pour une option Basket (en dimension $d \\ge 2$) on approche le prix par Monte Carlo mais on peut utiliser des approximations pour construire un problème unidimensionnel proche du produit Basket. Ces approximations servent de variables de contrôles: **on ne rajoute pas une erreur, on retire de la variance**.\n", "\n", "On rappelle que, en posant $\\mu_i = r - \\frac{1}{2}\\sigma_i^2$,\n", "\\begin{equation*}\n", " X = \\biggl(\\frac{1}{d} \\sum_{i=1}^d S^i_0 e^{\\mu_i T + \\sigma_i \\sqrt{T} \\tilde G_i} - K\\biggr)_+\n", "\\end{equation*}\n", "et en introduisant $a^i_0 = \\frac{S^i_0}{\\sum_{j=1}^d S^j_0}$ (t.q. $\\sum a^i_0 = 1$) et $\\bar S_0 = \\frac{1}{d} \\sum_{i=1}^d S^i_0$ on a \n", "\\begin{equation*}\n", " X = \\biggl(\\bar S_0 \\sum_{i=1}^d a^i_0 e^{\\mu_i T + \\sigma_i \\sqrt{T} \\tilde G_i} - K\\biggr)_+.\n", "\\end{equation*}\n", "La variable de contrôle proposée est obtenue en échangeant l'exponentielle et la moyenne pondérée par les poids $\\big(a^i_0\\big)_{i=1,\\dots,d}$:\n", "\\begin{equation*}\n", " Y = \\bigl(\\bar S_0 e^Z - K\\bigr)_+\n", " \\quad \\text{avec} \\quad \n", " Z = \\sum_{i=1}^d a^i_0 \\big(\\mu_i T + \\sigma_i \\sqrt{T} \\tilde G_i\\big) \n", "\\end{equation*}\n", "La variable aléatoire $Z$ suit une loi gaussienne $Z \\sim \\mathcal{N}(m T, s^2 T)$ avec\n", "\\begin{equation*}\n", " m = \\sum_{i=1}^d a^i_0 \\mu_i\n", " \\quad \\text{et} \\quad\n", " s^2 = \\sum_{j=1}^d \\Big( \\sum_{i=1}^d a^i_0 \\sigma_i L_{ij} \\Big)^2. \n", "\\end{equation*}\n", "Ainsi l'espérance de la variable de contrôle $Y$ est connue par la formule de Black-Scholes, car elle correspond au prix d'un call de strike $K$ d'un actif Black-Scholes de dimension 1, de valeur initiale $\\bar S_0$, de taux $\\rho = m+\\frac{1}{2} s^2$ et de volatilité $s$ (à un facteur d'actualisation près... attention à ça). On a donc \n", "\\begin{equation*}\n", " e^{-\\rho T} \\mathbf{E} \\big[ Y \\big] = P_{\\text{BS}}\\big(\\bar S_0, \\rho, s, T, K\\big),\n", "\\end{equation*}\n", "où \n", "\\begin{equation*}\n", " P_{\\text{BS}}\\big(x, r, \\sigma, T, K\\big) = x F_{\\mathcal{N}(0,1)}(d_1) - K e^{-r T} F_{\\mathcal{N}(0,1)}(d_2),\n", "\\end{equation*}\n", "avec $F_{\\mathcal{N}(0,1)}$ est la fonction de répartition de la loi normale centrée réduite et la notation \n", "\\begin{equation*}\n", " d_1 = \\frac{1}{\\sigma \\sqrt{T}} \\Big( \\log\\big( \\frac{x}{K} \\big) \n", " + \\big(r + \\frac{\\sigma^2}{2}\\big) T \\Big)\n", " \\quad \\text{et} \\quad\n", " d_2 = d_1 - \\sigma \\sqrt{T}\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "id": "47e5e558-f825-406d-8175-7f6cb8c427e0", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: préliminaires pour la variable de contrôle\n", "\n", "- Définir la fonction `price_call_BS` qui code la fonction $P_{\\text{BS}}\\big(x, r, \\sigma, T, K\\big)$ définie ci-dessus.\n", "- Initialiser les paramètres $\\bar S_0$, $(a^i_0)_{i=1,\\dots,d}$, $m$, $s^2$ et $\\rho$.\n", "- Calculer $\\mathbf{E}[Y]$ par la formule fermée.\n", "- Calculer $\\mathbf{E}[Y]$ par un estimateur Monte Carlo à partir de réalisations de $( G_1^{(j)}, \\dots, G_d^{(j)})$, $j \\in \\{1, \\dots, n\\}$.\n", "- Vérifier que tout est cohérent." ] }, { "cell_type": "code", "execution_count": null, "id": "c91214d2-5864-4020-95dd-2bdee8445d3c", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 21, "id": "5d6c6dd5-4501-40b3-989c-0417d3912b36", "metadata": {}, "outputs": [], "source": [ "K = 100" ] }, { "cell_type": "markdown", "id": "60e52cfa-01e9-451c-93e0-9a3a4d900485", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: MC avec variable de contrôle\n", "\n", "Implémenter l'estimateur de Monte Carlo avec variable de contrôle pour le calcul de $\\mathbf{E}[X]$ c'est à dire \n", "\\begin{equation*}\n", " \\mathbf{E}\\big[ X \\big] = \\mathbf{E} \\big[\\psi(G_1,\\dots,G_d,K) - (Y - \\mathbf{E}[Y]) \\big],\n", "\\end{equation*}\n", "où $Y$ est la variable de contrôle introduite précédemment et $\\mathbf{E}[Y]$ est calculée par la formule fermée. \n", "Comparer les ratios de variance pour les différentes valeurs de $K$." ] }, { "cell_type": "code", "execution_count": null, "id": "f154c42e-2c8b-402a-b7fd-7dc9e8828fc8", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.5" }, "toc-autonumbering": true }, "nbformat": 4, "nbformat_minor": 5 }