{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Training \"inspecteur modèles\" BNP Paribas - Monte-Carlo methods - December 12-17/2019" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# TP 2 - Variance reduction methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1. Control variate\n", "\n", "We consider a stochastic model defined by $Y = f(X)$. The goal is to estimate $m = \\mathbb{E}\\bigl[g(Y)\\bigl]$ for a given function $g$, such that $g(Y)$ is a square integrable random variable.\n", "\n", "We are given a simplified model $Y_c = f_c(X)$ such that the value of $m_c = \\mathbb{E}\\bigl[g(Y_c)\\bigl]$ is explicitly known and such that $\\mathbb{E}\\bigl[g(Y_c)^2\\bigl]<\\infty$.\n", "In practice, $f_c$ can be a function easier to evaluate (with lower computational cost) than the function $f$, or simply a function such that the expectation $m_c$ can be evaluated explicitly.\n", "\n", "We denote $(X_i)_{1\\leq i\\leq n}$ a sequence of i.i.d. copies of the input random variable $X$, and we set\n", "\n", " \\begin{eqnarray*}\n", " I_n&=&\\frac{1}{n}\\sum_{i=1}^n g(f(X_i)),\\qquad\n", " I_n^c = m_c + \\frac{1}{n} \\sum_{i=1}^n \\bigl( g(f(X_i)) - g(f_c(X_i)) \\bigr) .\n", " \\end{eqnarray*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 1 (theoretical):\n", "Check that $I_n$ et $I_n^c$ are unbiased estimators of $\\mathbb{E}\\bigl[g(Y)\\bigl]$, and compute their variances." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 2 (a):\n", "\n", "We assume that the input random variable $X$ follows a uniform distribution on $[0,1]$, and that $f(x)=e^x$, $g(y)=y$, and we choose $f_c(x)=1+x$.\n", "\n", "Simulate the two estimators $I_n$ and $I^c_n$ and their $95\\%$ confidence intervals.\n", "\n", "Plot the trajectories $n \\mapsto I_n, I_n^c$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sample size = 1000\n", "True expectation E[g(Y)] = 1.718 \n", "\n", "Standard MC estimator: 1.714\n", "Confidence interval(95%) = [1.684, 1.744] \n", "\n", "Control variate estimator: 1.715 \n", "Confidence interval(95%) = [1.702, 1.728] \n", "\n", "Variance ratio = 5.6 \n", "\n" ] } ], "source": [ "N = 1000 # Sample size\n", "\n", "integers1toN = np.arange(1,N+1) \n", "\n", "############################################\n", "# Compute the exact value of m_c:\n", "# here, m_c = 1+ E[X] with X uniform on [0,1]\n", "m_c = 3/2\n", "############################################\n", "\n", "##################################################\n", "# Draw the samples of the random variables Y = f(X)\n", "# and Y_c = f_c(X)\n", "##################################################\n", "X = np.random.uniform(low=0, high=1, size=N)\n", "\n", "Y = np.exp(X)\n", "\n", "Y_c = 1 + X\n", "############################################\n", "\n", "############################################\n", "## Evaluate the two estimators, their empirical\n", "## variances, and their confidence intervals\n", "\n", "I = np.mean(Y)\n", "emp_variance_MC = np.var(Y)\n", "\n", "I_c = m_c + np.mean(Y - Y_c)\n", "emp_variance_control = np.var(Y - Y_c)\n", "\n", "radius_confidence_int = 1.96 * np.sqrt(emp_variance_MC / N)\n", "radius_confidence_int_control = 1.96 * np.sqrt(emp_variance_control / N)\n", "\n", "############################################\n", "# The gain in terms of variance\n", "variance_ratio = emp_variance_MC / emp_variance_control\n", "############################################\n", "\n", "############################################\n", "# The true value of E[g(Y)] for comparison\n", "############################################\n", "Exp_gY = np.exp(1.) - 1.\n", "\n", "############################################\n", "# Plotting\n", "print(\"Sample size = %d\" %N)\n", "print(\"True expectation E[g(Y)] = %1.3f \\n\" %Exp_gY)\n", "\n", "print(\"Standard MC estimator: %1.3f\" %I)\n", "print(\"Confidence interval(95%%) = [%1.3f, %1.3f] \\n\" \\\n", " %(I - radius_confidence_int, I + radius_confidence_int))\n", "\n", "print(\"Control variate estimator: %1.3f \" %I_c)\n", "print(\"Confidence interval(95%%) = [%1.3f, %1.3f] \\n\" \\\n", " %(I_c - radius_confidence_int_control, I_c + radius_confidence_int_control))\n", "\n", "print(\"Variance ratio = %1.1f \\n\" %variance_ratio)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXecVNX1wL/3TdmZ7Z1ddpe29N5BQQUbikQUQyIqxqg/LInExBK7RhPFJBpLsJKYmBh7wdhiAUVFpEjvbYGFhYXtbWan3N8fZ2dn+y7NRbjfz2c+M++9+967r8w595x77rlKa43BYDAYTjys9q6AwWAwGNoHowAMBoPhBMUoAIPBYDhBMQrAYDAYTlCMAjAYDIYTFKMADAaD4QSlVQWglMpSSs1XSq1XSq1VSv2qiTKXKqVW1XwWKqUGHZ3qGgwGg+FIoVobB6CUSgfStdbfKaVigGXABVrrdXXKnAys11oXKaXOBe7TWo86mhU3GAwGw+Fhb62A1joPyKv5XaaUWg9kAOvqlFlYZ5dFQOYRrqfBYDAYjjCtKoC6KKW6AEOAb1sodhXwYTP7zwBmAOBmmKI7Q/vGHUwVDAaD4YRm2bJlB7TWKUfiWK26gGoLKhUNfAH8QWv9VjNlxgNPAWO11gUtHq+j0hEpc/GsPP8gq2wwGAwnLkqpZVrr4UfiWG2yAJRSDuBN4KUWhP9AYA5wbmvCP7xTsI3VNBgMBsORpi1RQAr4G9LJ+2gzZToBbwHTtdab2npyrQJtLWowGAyGI0xbLIAxwHRgtVJqRc26O4BOAFrrZ4B7gCTgKdEX+NtkohgLwGAwGNqNtkQBfQWoVspcDVx90Gc3FoDBcMj4fD5yc3PxeDztXRXDUcDlcpGZmYnD4Thq5zioKKAjjrEADIZDJjc3l5iYGLp06UKN5W04TtBaU1BQQG5uLl27dj1q52nfVBDGAjAYDhmPx0NSUpIR/schSimSkpKOunXXrgpAGwvAYDgsjPA/fvk+nq2xAAwGg+EEpZ0tAKMADIYfMkoppk+fXrvs9/tJSUlh0qRJtes+/PBDhg8fTp8+fejduzc333xze1TV0ATtbAEYF5DB8EMmKiqKNWvWUFVVBcAnn3xCRkZG7fY1a9bwy1/+kn//+9+sX7+eNWvW0K1bt/aqrqEBxgVkMBgOi3PPPZf3338fgJdffplp06bVbvvjH//InXfeSe/evQGw2+1cf/317VJPQ2NMGKjBcBxw442wYkXr5Q6GwYPhscdaL3fxxRdz//33M2nSJFatWsWVV17Jl19+CYgFcNNNNx3ZihmOGKYPwGAwHBYDBw4kJyeHl19+mYkTJ7Z3dQwHgbEADIbjgLa01I8m559/PjfffDOff/45BQXhXJD9+vVj2bJlDBpkJgk8FjF9AAaD4bC58sorueeeexgwYEC99bfccgsPPvggmzZJjshgMMijjzaZU9LQDrSrBaDjt7Xn6Q0GwxEiMzOTX/2q0XThDBw4kMcee4xp06ZRWVmJUorzzjuvHWpoaIr2dQE5Ktv19AaD4fAoLy9vtG7cuHGMGzeudnnSpEn1xgUYjh3a1wWE6QMwGAyG9qJ9FYBl+gAMBoOhvTBhoAaDwXCCYlJBGAwGwwmKcQEZDAbDCUpbJoXPUkrNV0qtV0qtVUo1ivVSSvVWSn2jlPIqpQ4i1Z+xAAwGg6G9aIsF4Adu0lr3AUYDv1BK9W1QphCYCfz5YE6uMRaAwfBDZ+/evVx88cVkZ2fTt29fJk6cWDvw62B47LHHqKw8+NDwcePGsXTp0oPez9AGBaC1ztNaf1fzuwxYD2Q0KJOvtV4C+A7q7EofVHGDwXBsobXmwgsvZNy4cWzdupV169bx4IMPsm/fvoM+VksKIBAwjcWjwUH1ASilugBDgG+PyNlNFJDB8INm/vz5OBwOrr322tp1gwcPZuzYsdxyyy3079+fAQMG8OqrrwLw+eefM27cOH784x/Tu3dvLr30UrTWPPHEE+zZs4fx48czfvx4AKKjo7nnnnsYNWoU33zzDZ999hlDhgxhwIABXHnllXi93na55uOJNo8EVkpFA28CN2qtSw/lZEqpGcAMANIBFURrMNOaGgyHx40f3ciKvUc2H/TgtME8dk7LWebWrFnDsGHDGq1/6623WLFiBStXruTAgQOMGDGCU089FYDly5ezdu1aOnbsyJgxY/j666+ZOXMmjz76KPPnzyc5ORmAiooK+vfvz/3334/H46FHjx589tln9OzZk8svv5ynn36aG2+88Yhe84lGmywApZQDEf4vaa3fOtSTaa2f01oP11oPl7MH8PsP9WgGg+FY5auvvmLatGnYbDY6dOjAaaedxpIlSwAYOXIkmZmZWJbF4MGDycnJafIYNpuNiy66CICNGzfStWtXevbsCcDPfvYzFixY8L1cy/FMqxaAkqnp/was11of2TR+KojPBw7HET2qwXDC0VpL/WjRr18/3njjjUbrtW6+fy8iIqL2t81mw99MK9DlcmGz2Vo9nuHQaYsFMAaYDpyulFpR85molLpWKXUtgFIqTSmVC/wGuEsplauUim31yCpoLACD4QfM6aefjtfr5fnnn69dt2TJEhISEnj11VcJBALs37+fBQsWMHLkyBaPFRMTQ1lZWZPbevfuTU5ODlu2bAHgX//6F6eddtqRu5ATlFYtAK31V0CLXnqt9V4g86DPXmMBGAyGHyZKKd5++21uvPFGZs2ahcvlokuXLjz22GOUl5czaNAglFL88Y9/JC0tjQ0bNjR7rBkzZnDuueeSnp7O/Pnz621zuVy88MILTJ06Fb/fz4gRI+p1PBsODdVeppXqqDTnjGHvQ1/RoUO7VMFg+EGzfv16+vTp097VMBxFmnrGSqlltf2oh0k75wLSVFe3aw0MBoPhhKXdFcA//9muNTAYDIYTlnbPBrp5c7vWwGAwGE5Y2nlGMBPaZTAYDO1Fu1sABoPBYGgf2r0PwGAwGAztQ5tzAR0djAVgMPxQKSgo4IwzzgAkJbTNZiMlJQWAxYsX43Q627N6h828efOIjIxk9OjRB1Vu9uzZxMfHc+mll34f1Tws2lcBGAvAYPjBkpSUxIoVkoDuvvvuIzo6mptvrj8flNYarTWW1c7djYfAvHnzSE5ObpMCqFvuF7/4xfdRvSNCu/cBaA1z5sDnn7drTQwGwxFiy5Yt9O/fn2uvvZahQ4eya9cu4uPja7e/8sorXH311QDs27ePKVOmMHz4cEaOHMmiRYsaHc/v9/Ob3/yGkSNHMnDgQObMmQPA66+/zoQJEwDYvXs3PXv2JD8/nzlz5nDhhRcyYcIEevXqxe9///vaY/3zn/9k5MiRDB48mOuvv55gULwQ77//PkOHDmXQoEGcffbZbN26lTlz5vCnP/2JwYMHs3DhQubOncuoUaMYMmQIZ599Nvn5+U2Wu+uuu3jsMcnN9N133zFq1CgGDhzIRRddRElJCQBjx47ltttuY+TIkfTq1YuFCxcehSfROu3sAtIoBf/3fzVLxiAwGA6NG2+EFUc2HTSDB8Njh5Zkbt26dbzwwgs888wzzSZ7A5g5cya33noro0ePJicnh0mTJrFmzZp6ZZ577jlSU1NZvHgxXq+X0aNHc/bZZzN16lTefPNNnnnmGebOncsf/vAHUlNTAXFBrVmzBqfTyYgRI5g0aRJ2u523336bhQsXYrfbmTFjBq+88gqnn3461113HV9++SWdO3emsLCQxMRErr76apKTk2tTThcVFXH++eejlOKZZ57hkUce4eGHH25U7oMPPqit+2WXXcZzzz3H2LFjueOOO3jggQf4859l4kStNYsXL+bdd9/l/vvv56OPPjqke304tLsLyAh9g+H4Izs7mxEjRrRa7tNPP2Xjxo21y0VFRVRVVeF2u2vXffzxx6xfv55XXnkFgJKSEjZv3kynTp2YPXs2/fv359RTT2Xq1Km1+0yYMIGEhAQALrjgAr766iv8fj9Llixh+HDJolBVVUVWVhZut5vx48fTuXNnABITE5us686dO/nJT37C3r178Xq9tampm6OgoACPx8PYsWMBSWE9ffr02u1TpkwBYNiwYc2mxD7aGAVgMBwPHGJL/WgRFRVV+9uyrHrpnD0eT+3vUCu4pQ5jrTVPPfVUbYdzXXJzc7HZbOzduxetNapmdinVYJYppRRaa6688koeeOCBetveeuutRuWb4he/+AV33HEHEydO5NNPP2XWrFktlm8tz1ooLXZLKbGPNu3cMxPETPVpMBzfWJZFQkICmzdvJhgM8vbbb9duO/PMM5k9e3bt8oom3FgTJkzgqaeeqhWSGzdupKqqCp/Px5VXXslrr71Gt27dePzxx2v3+fjjjykuLqayspK5c+cyZswYzjzzTF577TUOHDgASAt9586djBkzhnnz5rFjxw4ACgsLgcbpqUtKSsjIyEBrzT/r5LBpLo11cnIybre71r9/LKawbncLwMwHYDAc/zz88MOcc845dOrUib59+9bO5zt79myuu+46XnjhBfx+P+PHj6+nEACuueYadu7cyeDBgwFITU1l7ty5PPTQQ5xxxhmcfPLJ9OvXj5EjRzJx4kRAOlkvueQStm7dyvTp02v3vffeeznzzDMJBoM4HA6eeeYZRowYwdNPP83kyZPRWtOxY0c+/PBDJk+ezNSpU3nrrbeYPXs29913HxdeeCGZmZmMHDmSvLw8gEbl6vKvf/2L6667jqqqKrp3784LL7xwVO/zwdK+6aAv6sUFuRt45x1ZZ9xBBkPbMemgm2bOnDmsWbOmNhLnh8zxnQ4aMyGMwWAwtBft7gKqsQQNBoPhiBAaY2BonXbPBmosAIPh0DGTpR+/fB/PtlUFoJTKUkrNV0qtV0qtVUr9qokySin1hFJqi1JqlVJqaJvOroJHfOyKwXCi4HK5KCgoMErgOERrTUFBAS6X66iepy0uID9wk9b6O6VUDLBMKfWJ1npdnTLnAj1qPqOAp2u+W0ZBzchog8FwkGRmZpKbm8v+/fvbuyqGo4DL5SIzM/OonqNVBaC1zgPyan6XKaXWAxlAXQUwGXhRS1NkkVIqXimVXrNv85j5AAyGQ8bhcNC1a9f2robhB8xB9QEopboAQ4BvG2zKAHbVWc6tWddw/xlKqaVKqaWyxpiuBoPB0F60WQEopaKBN4EbtdalDTc3sUsj6a61fk5rPbw2htWkgzYYDIZ2o01hoEopByL8X9Jav9VEkVwgq85yJrCn9SOHFUDDVBx+v/QPBAIQEQG7d0N0NHTq1JYaGwwGg6E12hIFpIC/Aeu11o82U+xd4PKaaKDRQEmr/n+o1wdgs9Xf9LvfQXIydOgA2dnQrx/UJOszGAwGwxGgLRbAGGA6sFopFQravAPoBKC1fgb4AJgIbAEqgZ8fbEXqKgC/H+rM4UBBwcEezWAwGAyt0ZYooK9o2sdft4wGDn4etDoWQLBOQNBbTTmZDAaDwXBEOWYm6vT5wsngysvbty4Gg8FwItDucwLXJZQa2qSHMBgMhqNPu+cCqktI8NfMx9Ak27cfxeoYDAbDCUQ7WwD1F0MKYN++5neZMePoVcdgMBhOJNo5HXR9F9Du3dC3L8THyweguLj+LpWV31PdDAaD4TjnmHIBzZgBe/bAunUwaBBUV9cvPW0a7N37PVbPYDAYjmPaUQEoGiqAr78O/05Ph1tvrb9HZibk5tYPGTUYDAbDodHOfQDN5wJKT4d774W8OuOJs7LEKjDZbw0Gg+HwOSbGAfTv33hdWpp8p6aG12XVZBvKzT36dTIYDIbjnWPCAlizpvGm9HT5turUMKQAduw4yvUyGAyGE4B2tgCad+YnJoZ/L10KmzdLhJDbDZ9/fvRrZjAYDMc77agANNj8zW6Njg7/HjYMuncX4T9mDLz+ejhthMFgMBgOjWMqDLQuUVFNrx82TEJB3377KFXJYDAYThB+cAogFBr62WdHoToGg8FwAnFMRAE1RXMKIDERRoyAbdu+3/oYDAbD8cYPzgIASEiAoqKjUB2DwWA4gTimksHVpSUFEB9vFIDBYDAcLsesBRAR0fxexgIwGAyGw6ctk8L/XSmVr5RqYrgWKKUSlFJvK6VWKaUWK6WaGNfb3MEbK4CHHoJFi0C1YB2EFIAJBTUYDIZDpy0WwD+Ac1rYfgewQms9ELgceLztp28swTt2hFGjWt4rJUVmD2uYKtpgMBgMbadVBaC1XgC0MEcXfYHPaspuALoopTq06ey2ALjrH9pma323zEz5NjmBDAaD4dA5En0AK4EpAEqpkUBnILOpgkqpGUqppUqppbUrT3qkXpm6I4CbI5QTaNeuQ6qvwWAwGDgyCmAWkKCUWgHcACwHmszxoLV+Tms9XGs9vHalz12vTFsUQI8e8r169aFV2GAwGAxHQAForUu11j/XWg9G+gBSgLZP3e6Nq7fYUvRPiORkUQJz55qOYIPBYDhUDlsBKKXilVLOmsWrgQVa69I2HyBQf1rits756/PBN9/Ac8813vbccxJFVFXV5loYDAbDCUdbwkBfBr4BeimlcpVSVymlrlVKXVtTpA+wVim1ATgX+NVB1cBRX0pXVLRtt9BMYU2lhr79dvk2fQQGg8HQPPbWCmitp7Wy/RugxyHXoIECaDgRfHP83//BX/8qlsDbb8OFF8r6tWuhsCawaNcu6NnzkGtmMBgMxzXtnwzOfmgK4NFHpb/gzTdhyhRRAFu31p9e0lgABoPB0DztrwAcnnqLbVUADkf9DuB33pFJY+piFIDBYDA0T/srgEO0AECUQEsYBWAwGAzN0/4KwOaBoc9Dl8+Bg1MALREdbRSAwWAwtET7KwC7B86fAVeMBw5OAZx/vnzfd1/99ffeC2eeKQpg4ULweo9MVQ0Gg+F4ov0VQOraeosHowD+/nfYvFkEfkxMeP1vfyvfa9fKJPI33HAE6mkwGAzHGe2vANJX1ls8GAXgcoU7fkNJ5HbsALcbunQJl3v+eSgrg5wcmDULgsHDqrHBYDAcF7S/AmjAofYBPP00pKVBerosP/RQ/e2xsdC1qwwSu+UWsQ4MBoPhROaYUgCWdegK4OKLZXRwKDLI5ZKJY5ri0Udh+vRDO4/BYDAcLxxTCsDhaLnDdt482LbtyJxr+XLo16/tqScMBoPheOOYUgDx8c1bAFu3whlnwOWXt+1YwaD4/QGcNanqfvc7OOuscJl16+DJJ+X3Z59BdjasX39odTcYDIYfGseOAtCKoNbNKoA31syFrIV8/XXbDrdrl0wbCWG30OWXw8cfy/pQwrjbb4d9+6RfYNs2eOGFw7sMg8Fg+KHQajK47w2l2V/oxeNxAXDRPa+zIfpp+nVNpntidx5a8RBcBdzXtgkAdu6U706d5Pezz4Yjg2w2+MMf4PHHJf10Wlp4vy++kARzrY0yNhgMhh86x44FAPCz8ezYAWvz1/KW7Sesq5rP6+te56Gv6oT0xOxpNYzT44EPPpDfv/+9fKek1C+jFPz5z/XXPfYYLF5sxg0YDIYTg2NLAWQtYtEi+GTrvObLjL+bgoKmN4WSw11zjcT7A5x7rrT4ly0Ll3vxRekEvv76+vv/6leSZvof/4C9ew/5KgwGg+EHwbGhAOb/rt7i6t1bmi63bwB0+bx2Mpi6PPOMhJGWloqfH2R0cHIy9O0bVgCbNsHPfgZDh4b3jYmBP/1Jft94o7iA/vrXw7wmg8FgOMY5NhRASVb496+68tL6JuZ5BFhxBSRuY+2Oxs3zZ5+V788/D1sCHTrI97Bh8O23sn7z5vr7XXyxhJ5ed50s9+0Lo0bBp58e8tUYDAbDD4JjQwGoOk79hBy8QQ98cTe89hop9q7hbTvHAPDFzsYuopCPf9MmCATkd0iojxsHRUViDUyaVH+/Sy+V0NPExPC6M86AJUvEmmgrgYAknmuK/HxRMkcq06nBYDAcCdoyJ/DflVL5Sqk1zWyPU0r9Vym1Uim1Vin184OuRZf5jdftGQbrprL/rm1Q0AOWXwF7hkNxZz4qmI3W9aOBSkrk+4sv4MABCe/89a9l3eWXS4qI0FSRIH0Db78NQ4bIcnU1zJ8vQn/gQBlHsHx52y/hyScl8dx//wvjx8Mrr8h6n08sEZdLZjD76U9l3cqVNNuXYTAYDN8HbbEA/gGc08L2XwDrtNaDgHHAI0op50HVYtBL8l1Xppd1DP9+chPMfQG0jYiVv2SHXsg/Vvyj3iFCA7jee0++p02TSB+Q7//9L1z2mWckSuiCCyAjA26+WdaffrpYCdfWTHe/ZIl8r14Nt90WHljWFKHW//nnixtq2jQR8JdcUr/ca6/JwLTBg0UJGQwGQ3vRqgLQWi8AClsqAsQopRQQXVPW3+YaBOoMRVB11ld0aLK4fclNRHt78J81/6ldl5vbWDj37Vt/2V+nRjk59beFQkVBWuchS+HBB8WauOEGePhh6Txuju++a7wuORneeEN+//e/8MQT9bcvWSKD0AwGg6E9OBJ9AH8F+gB7gNXAr7TWTUbqK6VmKKWWKqWWhlcGmj5qRWqTqysrFNH7z+Db3G8JBGXff/9btkVGhsv95S/h3zfcUD/qZ9YsKC4OL0dENJ0iuqhI+ha++EKW330X9u8Xi6Auu3dLqopbbpFtnvrTHDN1qvQ93HCDKKurroKPPpJtZuSxwWBoL46EApgArAA6AoOBvyqlYpsqqLV+Tms9XGs9vM7a+oX2DIMl14Hf1ewJgztHU1ZdxvLcDTz/fDitw/jxNQVS1vHUptsIBIM8/HDTIZ0NO2yVgqVLZf3+/fC3v9Xffv/90tF7++3SR/DOO+FtIQVx8cXQv78olBdegFtvlcij114Ll83IgDlzYMIE6Zx+7jkzP4HBYGgfVMPO1CYLKdUFeE9r3b+Jbe8Ds7TWX9YszwNu01ovbvGYHZXmGiBIfTW07Gr47/MtVyhhK/yqO3z8R1h4S+3qefPg9DOCcO1g6LCauzO+5IH/G1tv12XLJCx06FB4+WXo2bPpU2zZAj16hJcLC6FbN+ls1lr227hRfo8dK/MLFBSEJ6ZpC6++KkrDZpPIo0sukb6B3r3Bbj+4YxkMhhMDpdSy+o3oQ+dIWAA7gTMAlFIdgF5A25M2N6xBVTiJ/513NrNPUTZUJMPZt0L3DyFWZn8fOhSGPD0cOoiP5oHX367d5Y03JMIn5Ar67jvo1Ut8/k2RnV2/H2H5ckhKCo8x2LRJUkns3i1Ww89/fvAC+8IL5TsQkMFrV1whCsDlEkth4kR4//2DO6bBYDC0lbaEgb4MfAP0UkrlKqWuUkpdq5SqiZXhAeBkpdRq4DPgt1rrA4dUGw2osEUycmQLZddNle/LJsJvOhF74ymc+8bJLN9bJ3bz5EchTZbPOis8b3BcXLiI0wl79jQ+vFLSqt+5U36/9BJERYWF/KhR4vMPKamf/vTgLjV07pB7atYscRmF2LcPPvxQ+g7W1ATgPvus1OWCC5pXXAaDwdBW2uQCOionDrmAGrJtPLwoA73+85/GYZS12Lxwd9P9BGd1O4slf/0VxRMnEbntYhbd/DIDBoS3r1snmUCfqzPgWGsZsBUIhKeVDHHrreFUEdddJ9NPNqSyUuYiPlxKSkQpvPQSjBgheYsATj65fr9FdrbUZdgwUSRZWfIxGAzHN8eaC+jIkliTByj7Y77avLLJIuPHA4EIWHCHrNgcDqgfmj6UDy79gC0fnId9+bVUdnuFzfa36u3ft284dUQIpWTAVu/ejc93223SsRsZCTNnNl3twxH+Xq/U57rrxDV1xx2ipP75z/CAspDwX7AAUlMl6ujmm+VejBkjiuCee2SWswkTwqOhDQaDoTmOPQsA4H4v3BMhv+vm/09ZB4XZIvyBjhlBRo+t4q1Xo7jsL88zt+om8m/Jx2UXyyC3NJesv2QRGxHL1plbSY5MrneanBwZIBYa+FV7+vsl5r9Tp/C6oiKZsUwpifo577z600lmZEgiuT17xB00alT9Y+7aBZMny3E6dYIzz4S775Zt06eHQ1lBwkTtdrE6Qm6rJ5+EAQMkcqi0VBREQYG4iJQK902EuOkmiI0VxfLkk5LhdORIOXZDC+eHTMjyUqr1sgArVkiU15lntn0fg+FY4khaAGit2+VDOpr7mvl0WhD+jZZPwhZZvvSc2nU5OVrPnCm/X3xR16O4WOun315ee5zr37teD3x6oO43u59+YtETOrckVz+44EF997y7da+he8PnqflERekW8fm0LinR+ne/0432Ba2vvlrrFSu0Dga13rpV6/T0xmWeeEK2gdZKad23b9PHOuUUrf/wB60rKrT2++vXY9as+mU7dZJjNXWcpj7nnad1eXnL13qssXGj1jNmaH3OOXKt6elaP/CA1h6P1nfdpfWzz8rzqazU+oMPtL79dq0vvFDr3r21ttnkui+4QOtHH9W6oECekcHwQwFYqo+QHD42LYC3/w4XXinl/rYIm47A794Dl54n2+8LAoqVK+Gyy2Tw1VUPfsKOzD8y9+K5RDoiueYaeC7nZjj5kVbrMjpxIpMq3qLf2Yu5cNhYQkOSV6+WuP7WmD0bfvnL1suBDAYLzUNcl+nT4V//an3/sWPFQpk3T+YvGDtWIpHKy9t2/qaIjpYO57E1EbNaS/+L0ynXlpQkaTIuu0ysn1NOgYSElo95OBQUiLVlWTIvQ16euOZuvVUG4+3aFS7bs6eE7LZ1LIXdLtcSShluWbIulBCwZ0+5D/fdJ53+rREImHBdw/fLkbQAjk0F8O31MOqp5nd+fCsUdcPlqjPq9urRkPktT579LL88aQYTJsDH2QPI7OggN9B6VjcLiyBBruhyL+/fch/794u7ZMoUCcdsCp9PhPYjj8C6TVXQey6su4gBI4tZvTgBgjVpLjqshPJ0zp0YYN3idKqqpMO5Kf7yF+nMveUW2L5d1k2dKrmL6rqcINwx/PzzMq3l5MkivO32xlFCV1whbqCGREeHlcezz4pwHzmydYWydKn0OwSDIkS1PnSXyuLFcoyhQ8V9Fkqf0RwZGZLgz2aD7t1FMTR3P7t0ESGdlQWrVtW/ruRkSfURoq4rzeWSkNxf/EK+p02TaKyLL5Y+lj//GTp2hM8+kzr5w0lUAAAgAElEQVT//e+yz9GiqKh1pVtWJgkNfT4JIIiMlIi30PSmDZ+R3y9BB1FRcv+dB5fBy9BOHP8KoDgL4nc1sxF4+x+wskFinl91hYQcADoFT2Pnc4/AtcMZV/UEP5rs4+55d1PprwTAqo4nuGcQOCogYykNiXJEMWDDq6x551zKyywmT4Y332zc0rvpJnj00ZqFs26FMX+C8g4QvQ8rdwwZJT+m9/jlfJL/Yv0dl/8cMhfBwpthuVg6V18trc6MjHCxlSulYzotTQT9mDFN346f/CQ82viTTyTktSn69xeBWXcUc2tMnixRTw891LTlEkIpEbL9+0OfPhJCu2CBCJk775QO6hAffyxK9aqrZHxFKOleXRyOsBJLSRG/fUSEdJjHxsogvbqzvL3xBgwaJNe3f7+cY+rUxkJt48bwvbrrLhGAW7eKZXPTTW2/Lw0ZPVr6cbKzZYxIZaUohIIC6cjfulXKrV8v1oZSEuo7YoSUGzdOFNWyZfLdqZP09SQnh0OFr7lGlFGXLtLX5PHA11+LEl+6tOl8VP37i0JfulT6jfx+yMwUoe9yyTsG4QZEUpK8TwMGiFIwHHsc/wpAK1Aa9dUd6LEP1t8WsIvwf3dOeF3sLvhNJyhPhegGTcEv7kCd9hC6YcqJr2+Bb37Nj6d5eCNhOHhjIGFH/TIbfwSf3wsRZTxw1Tjuugs8fg/vb3qf0QmTyexY08KP3A83dgVnBew8GTo1nhjAsekn+FQp9Pio/oaP/8RjP/0NF1+5n9SoVFQTzeiKCglZLS4W909z/PSnIty9Xlk+55xwzqGWaNiJfPrp0toNtapPO02snMMRkjabRDe5XE0P8IuPl+sbPRoWLZJ1XbuKEKybNvvii0WQbd8uIcJDhsBFF4mgPFz27xfBGAyK2+nNN+Gpp0RQZ2VJa7q4GLY1McwxdA/rWlTNkZUlLq2UFBHMmzeH739dxVcXp7P1+STOOAN+9COpw7vvimtw7Vo59ogRojTfe0+EvWWJVTFwoNz7HTtEEYcs6tRUURS9eomrcdMmUdZTpoiLbMsWeT6hRpHfL79Nx/rR5/hXAADVkfBgOUw/C7I/C6/XQEFP+MfnUJ4GKOyn/wH/qXfBPz+DiydDRPP/wMmVc5kbOTm84svb4NtfQkUaxOeQ3dNHhytmsnDvJ/Wzk85eww2/tLEo+jaWlM2FLRPg9degOgr7zAH44zcw6Kt1/OPPvWWOgT5vQe+3id9+JcVFNthxKgCLl/m4ffYiPvvPAJh0DfQPJwrKiOrM0PShFFXvZ3CHwThtTn7S72Jm3TCittVuWeJ6eeABEU7r14vw++abcFVnzpRIn1mzxFUxZYpEKL37rrR4m4oamjJFWq0ffSSupd/8RiyO3/42PBYhRMeO4cFzTqecZ+rUxuMQrr5aWpK33y7HDtGpk1gJdVN0Z2XJyOhQxtQ775QoqYgI8fn/5z+y349/LC6uYLD9fO8HDojwXLdO7te334rgDKUkj4sT4Z6eLlZJv34imFNTRZCWl0PnziJ0H35YrJWKCtnvd78TJXTGGbBhg0QtXXutXOt770k/xddfi8upqQGMIBZCKONtp05yH6+6Sp773r1St+bu3a5d4ka6/365zmCwcabd7GyxaEaMENfYunXyfNLSZIKls88O9ycZWiYQkAbPzp0SmZac3LIS9fnA6TwRFMDeAfDKXLixW4vHiV3wDP7RD1JZboMPH5c0EdMugKTw3I9xlUMoefFvkmG0LAM6L0Bdcj46omYWmaANDvSCoJ3o3RdSvvJ0uOgSiK35hyldM0q5zomDFijN2V0n8nHO+0zrcQ1/u/AZ3G75k116qXTU1uXKK+snmfvfJwHOefyXMOIZRkVeyrebNkNmEymUdo4hrWo8Kb23srrqfdg9ipNdVzPnDwOgoBd9+4Rt9d/9Tv6Qr74a3n36dGnFX301DB8uLb3168V9NaeOITV5Msyd2/j0ltW4k3XsWGmZh4TehAki0DZsECEVEyNC7fzzpeW4Y4cc50CDMeIulwi8uq18p1PGOHg8ci9Drp3ISGnFpqSIAPzsM5lT4dprw+6KwkJpAbfFn621tPqjo+Grr8S1tHmznGv06INvzX70kdRl9+5w+vFrrxVhnpAg2x95pL7ytdtFkNZV4CBjVdata/l8EyeKsHe7RXlu2yaCefFiOZ/TKcJ8xw55th6PKN1QmpE+feD660XJNiQYlE9lpbxLKSlyHc8+K8esrpY+lZCVOHy4HOfrr2W5Rw+5l2efLcpm2zaxJtLS5B3r0EGCCkaOlHtVXi7CLZRnKzNT7k15ubjF4uOlfwokh9ecOXLsmBix0HJypL6BgFhXp54qiqpjR9m/d295vnl54T4hkPMFAnKuENXV0lBKSZHtOTmi5PfvFyW9fbvc89JSeW/WrxcFe8YZcn+Li6XhU1Ym72PfvrL/woVinaemSt327JG5Q+qmqu/aVULM09Lk2B6P1M3hkGv46CNYv/44VQBjMk7j691fhFcsmgmjGyTRbyvVbhHcT60mxZHN/pHXiWvp/aeg/yvMvDOPlz5dRUHWP5ve3xODPSKIXzXoed07EP7zX0jaAj87AwCH5aD8jnKqfFWUV5eTESuO/LPPFp88yEvSs2djv+qePfX9/kSUQNZC0JZYOINeRI14Du2osWqCFgQdYPeG9/G5se8bRXZ0fzbu3QklnWHhTTz2xxhufGgl5IwDFETl89X/UhgzRlFVBc6IAI8/ZnHbbeqQU0vUddm0lcjI+hbB4TJkiMzpsGSJKMB+/UQQBQKicL77Tp5FQgL4g348fg/2YDS//GXjrK91GTFCWs8/+tHB1cfvl7rce2/4+YcYNEiUwtq14p8/6yxJMgginF56SZTo8uWS8qN3b2lIFBbKffvkE3mPQi3uttTlj3+UFr3PByedJMcKKe6+fSXwoEMHaX0uXiz9DF27tnxcCAvKuLiwwt2+XYIN3n1XriM5WVKgZ2SIYKyokOdTUiLrD4YOHeQehIIj2kLIddbQ6nW5RPEXFMj6uLjwZ98+UWxdu0rLvKVBldHRosTmz29sVdclJkb+K1FR0qjZtk1kwfjx8m6mpsozX7RIGjZVVVKXpCRRJBkZ8syysmDLluNUAQAQVGDV1GlfP+iwFnwR4PA2UbiNvPUvmDK9yU1DXRdRUFrOjqLdsH4KDJ0TbvmH0Ape+i+M+mvYh7/+AugjfhmbdhFQ4jy1BaM4qfRP/PriYdz58xFsWC/NyB1FO7li7hXklecxtc9Uxncdz2ldTsNSFvc/nsu9L35Ix1QXqZFprHhnHDabxcD+tpppKTV0nSeCvbB7WEkkboVecyGiFFLXgLO+VHXanFQHqsETJ/MuRJST7MjkgK/OP68iGTafh7XxQoLrzycyUuH1ystnt4fdAI3RuFyq3twH0dHyJwi1tFojMlKsk1dflRbqn/8s03SuWyetn6V1+udjYqQ15KMCXCXEDfqcYafuY799ORtXxVPt3A2lHaHbZ2AFYNdoiNkr98ZVBCqIZdcEo3PB8kL+APDGElU2BNeBkWR28tM1PYneneMpL4zhr3cPFCWM/ElHj5YWbHGxtEZTU+WPedYEP/Gxtib7bkBCif/0J1FSmZnSX9Fa5+rKlVImKkpcOKEw1W3bxIpcv176U37yExEQ0dGt3+vCQmmJx8XJMyotlRbpjBlNC+KUFCmTnS3uqtRU6XOJjpbzR0aG3XGhS1+/Xp5ldLQo+O3b5X4tXy7rYmPlGQ4ZItZBSNAWF8vzLSsTJZ2VJXXy++Xdi4mp/z45nVK/3btbv26nU45Rt5XdFpQK98eERKRSYuVERIgCi46W9yIlRe6T3S5W1aRJ8u4mJso937NH3Kkulxxr/Xo5Rkiwh+YXGTZM7msgIIoyMVEaEk6nWHZKiUL69a+PRwUQcrH4nWBv0NsVcIDNB9tPg65fcEh44qCoK6SvAGBY+jCW5S1rVEyVZaBjmnizNFCaDlEHwN6guVztBmdVM+eNJdJtp1I3nlTNrpz0SezH6vw1cn11zxWSJ0GFw5OBL78rAx0XUrXxNDYv7gpdvoDtp+M4+x58wx5vet+G9W/FpRHvHcQL0x7GQxGvrXuN7ond+d+axaz6ogtUdKBzeix9Ru3ko8WbIGMJ7DwFCrsTWzYKz/KLqK6omw9DE+ECr1cz7ZIAy1b42JcbWTt3cz2i9knlvDFw+p3QZT5xcTZK7Bsh6EB540mOd1OoNhHQDf7JQZsI/IOhxn1XN/FgQyxl4bRcRFRnULE/EX/sJlEmfheUZsh74ImDuF3EksXI2PP5xbnnMCyzH0v3LMVpc1LkKaLEU8Le8r0cqDxAdaAaT3WAb9fuJSkqjs5dNMmOTpTmpdAtPYnd26P5ZOF+dm2NAp8bLD8EHTjSNuGL2g75/aGkE3Z/PH6PGyJKiErfw1nneohM3UMxu9hTsk9Gjxd3oXOXAMnRcSS4E/D6vZR6S4l0RNZen6UsVu/Yw7x3U6kodcKB3sS54onydGffhmx5JkFbq2lFIiJEuDX5bFvA7RahaVnSGZ2VJa3uG24Q4e5wyLbOnSUA4dlnZfv8+WH3WESE9LUEg6KwSksbd6K7XGGl/dvfhqOnQhbMV19Jn86hYrPVtxIsK3xPiorkGhMTpZ45Oa3fp1DoblFRc9bH8agAakhzZ7K3qkGTpCJZ/nDAw0Pf5LeLLwW7p/HOAKVpRC/5AxWn3YC2V4InFlylsOF8+OAJUk99j/xev5fWYWssuVYEbdxOifCpg92XQCdXP7blltcqldawKWkp+oMtNEe0kpZnS0ItJMyDFljBNgn3o0JNtFZtnYI2+WDVKPEaIauVVLAyHryJshy/XbZbgZbrHno9Q2U8cRBwEBurqC5OwZPbAyqTIWccF121g9REF9ldXGwq3EjfDtl0TehKXlkeie5ENLBg9Vas2DyCBJg5aiZf7/yadfs3Ul5dytqcfeytyGNPfhXupAN4VQll1S1MBN2wni1cR4TlxtIOqvyVYsmpAKggxOwBWxuap9WRjSy8RgRtojgCEeAsA6WxbEGCBLCwiHRGUl5djkJhYREgIPfTUYWyAui6s/MF7FKvqgTYOwjcRRC9D/L7woHeEJ8DpVlikTrLZbk6Cg70EQV5oDcxyeX4ilPpnBHBvvJ9lBZGoj0x6Nhd4HPJffDW5Dqx/PLeKw1R+TXX4ZRjx+yR/3txV/DEQ/oyyPwW/BEy1sZWLfuH3KLVUVI3TxwU9AB3gWyPz4GkTVCVJJ9ABCRsIc6eTreY3gwaaOfUoR3o3dPOvG8O8N38zqxa76HQu5fYrFwKij2U5aWKN0IFpE7OCkjYJvehuIucy10IvkjQNql/VSJE1nR+hf7ffpdca0RpzaekpnGRJa5fX6RYrhFlco0+txwnZg9sPef4VQAdozPYU96gBd5SC7sJ7Hkn409vHIpZj6DFjcPu4lz3vUx48D4R8kNq+gPefxJ2nczg6a+w4nFJ2HPx7FkccC5mfs58AroJ4RxU4EmQh2mreZnLU+hrn8y6nPxad1Ety38G2oaj3/tEWJGUO7aTGZNJcVU55f5i7JYDp+WoHbvQKnUEUKw9ibL1J3H5aaewKvgKa/LX4AtKsyjKEUWFr6L54xwqR0gJKRROmwt/ZSQBW1nYGgzY5Rw2f/g82hLrsCwD+4FB+As7ipCqjoGdY3E6FQ6bnQp/Md06xjNxTBd+9qPufLF0P6PGVvK/rR/w/FdvsW+vHWJ2gd8N9ioo6gZViajyTHRRjZDbNwhKanoOo/dBVD7OxH1Epx6gsLyMqE5bqUidB/FbRSmHqI4Cb6zUP6IM3MUcMQI2LCKwVcfjK00QgeIsB2c5ygLtc4GrEBzNNJZAIurK0rDsAYLVEbJ/wAnOMpTND9qOjtojbrNAhGThjSgTYWXzgb1GeXkjRaE5WzhXU9RV8A2VfXPlT/RQ0/s4fhVAW4kkkUoKYdU0GPjywe2sFTy5EQp7iPvhljTIOU1e7pR1YjGE2HmSdLp2WdDoMLU+9jbSqXoCriW3synzTvkTpa0iLiKOEm/zNmF6dDp3n3Y3i3Yu5sU5UTDiqXAfSRPYLXuTFsa0/tNYuGshO0rqj3VQKO4ffz9oWLJnCUPSh1DmlVZvnCuOV9a8wvoD61s9fkNUTaJZzVGa7zJQE7Zh+UMnpOaE9ZdbozhTnoUVlIGBIddQa8fRSt4LT5y0vAMO8MZD7A5w17w/5YmisNx13qc2CDBbeRb4IolPqcLjrwJbNUH8WIEoKkpdUB2BWyVS7Q8QCPoheb28u831LVTFQHU0+KKlrs4ycJVIi9nmO3EE6vetPJr7mzZ8V0Pr2vruHg8KoFeE0s8eRlbKaFsi5YFCyQ6auFVaZ0GbmI8tmMoxEdGUeWsiampG7TZLA/9yTEQ0adHpJLoTsSkLm2Vn5Y7tlCIuK4flJNLpIsoRTZW/iqKqlh2LllIkuhOp9FWhlMJpc1DmLcMfrG9hKGw1glSeVXZiN+JdCbjtbsqryyn2FFHkKaJbQjdiI+Io9Zawt3wfSe5E4lzxOKxwjJtYAgqHZafK78Fld7X6vlX6KmUfmwOHZafCV0GVr4oKXyURtgj8QR8l3lKqfFUEtETZhHzMKZGpFHuKqfRVEuuKJT06HY/fQ0HlATTgsruwlKLK58Fm2XDanHj9Hip8lVhKEeWIotJX1UjRKiw0GhVwEuFw4AlWiFKwalwrh4BSCktZBOre/5CbS1vyjOx2tNb4gwGCTVmChiOKQjUexNke9Whi7Iy8L4pgKMGlalDfo1Tt8TtORAWglfjUrGZan1VJ4u8L2iHra1CaSJVEpS6gU0w3tFVNXvleFIokdyL7yvPrvVj9UvuRV1BKSlwcW4o2EdDNt+wjbG4Gpw9k494civ1hBTI68yRc9nAAeqXXx+JtG8QnWIsiwR3fqnKoJdQBHt699sWyasIvgjUvX4TNSWxELG67m9iIONwONw7LgdPmaNu5jiANLftA0I+9jiI6uGNpiqqKKawqIKBFuJd5SymvrnFlaaue0FdK0dR7bWkHURFuKrwegvgJNQltyo7dsuFyRODzB7HZNOXeKrSSd83CgQ4qtOWjxX913T6RmqsPP6P6SklhEWF31mxXeP1+AkFfnf3rNgkb7KlqnnmL9aD2+tCi2DSBBsJU/NGWskBb2JSNSGckaE25VxSy1tJHY7fZcDltRNgi0FpTHaimstqLy+HE7XBRHfDiC/oJ6iDVvgBKKXwBPxp/i53t9QjUvKe2Br24QTsOInHb3US5HcS6o2reJwc2y6q1SEMNp4rqCqoDPqq8fpS24a+pl6UUQeUnoH01/R0t3UNL7o86OCVvtxy4HS6cNic2y4aFRYQ9gqAO4PF7qPJ50Ghp7FgOLGXhC/oJaOmn8QcDBLSfYDCIRqMJorUmqJEMnmjG79A/fAVwSC6gL+6A0x5sfvuqS+GtfxOXXsgpV77Pe7OmEnNHL8psOwG4fvj1PLVUkszNGvI6ty36OeSM55yMSyjVuSyMvJ20zXfx3XMzmL3oec7tdSaFVQX0tp9Dz8sfh5FPNJmj6NTd7/D3Wf2Jd8XjC/pIiw4HaN94Izz+6aukTbsXZWnyqje1fI27R8DcOXDWb6HHR1hY9HVMonzZ+fQ/ZSuLN+ziwJIzCWa/R+KAxVQFS7EcPip8FdgtO267u1HH5amdT+Wc7HMo9hQzOG0w5/Y4l3hXfO32Sl8lFdUV7CrdRam3lOpANdkJ2WQnZgMivBbuWkhhVSELdiygX0o/RmWOwm7ZyYjJ4I11b/DautewKRv5Ffk4axRRr6RepEWnEemIZECHAeSV5fHv1f9mW9E2Kn2V2C07Se4kTul0CtMGTGNI2hDKq8tZumcpLyx/ga1FW8mMzSTJncSBqgN8tOWjNnXKhlphsY44ynxlaIKogBvtt8BRWdPx6WnZN34Q2P3x+O31fftWdTzB8iRwlmEpi2BZKpbDD44KgqoaKlKwpW4mYNXv27IFIgnYxIJV5enovQOhLF3cTZYPFZeHtjzS8Q3SKbpvINdNPIXzTs2gX5dU0tKaT0wXipDRuiZGPlBNfkU+NmWjQ3QHUQZHkPJy+PLrACu35vHRgv0s/84mo10tTZ/sSNxxFZRXl3HmyCzOnughOToehz+ef765j6pKG+X+IiZPiqB3VirRzmhyS3Nx290kRyZjs2yUeErYVboLS1lkxGQQGxFLQAewW3aqA9VU+arIKc7BZXcRYYug0l/J1sKtuOwuUh3Z7NruIjPdSaTLzu6SPN79ajPVVS7GnWpnY/kithVtp2RvAhtWR1NcZIEnnj3F+RC/A0dMCS6nRcCqolIXgKMK0BBRgTPSA1Y12PxUqzJAY9dughWJaJ8LR2QVOqIIH/L8xcoXZWNpJ7rajQq4AU2w2iV9ScoPWDB7w/enAJRSfwcmAfla60bJkZVStwCX1izagT5AitZNxD3W3a+OAlAVHdBRLbhiQpR2rI3Rd3oyqHbVdBZXJkFkAVZ1HOrRPQR//GN09w8BiCkfTFl0/Sgdh+XAF/RxafBjfnf5Wcz87618UPKnemWm9J7ChgMbWHdgHZcNuIzy9+/hHdcFkLJe/pANxwo0IC0qjUcmPMLk7J+QdMcIvAkr6oe4BiwoT4c3XoKsRbB0hvihvQmgLTp30dzy7IcMzOzBZT/KZGdepVg5jW5kADouheg86P0Otugizuh5Ej+aEI0n4MHj9/Dk4ifJrwjnSHLZXXSM6ci4zuOIsEfw0uqXKPWWNjr00PShDOwwkLkb5lLkadlisSs7GbEZVPmrcFgOXHYXW4u2NioX44whOTKZ/Ir8g+6MjnHGYAtE4wwkUEURwWo39kAMDrcH7SrEXp1EpS2PMl/jPhWbsoFWEkbakl8fwBOLrTIDe1k2BO1YkcVU6RKIKIKYvNYVhzcS/FFE+brhCZQTcOXjru6Eu7oTVnkW3gontv3DKF51kkRD9fhQomJs1bD1XCjuJAK/PJ1J51lMmQIp6V5W5+wmZ2VniisqSE5wkpgUpMewXfQcWIw34GHd/nW8s+EdVuevRmtNSlQKY7LGkBKVwt7yvWwv3k5uaS4V1RWUekuxlFXvuca74kmNSsWu7AxKG0SsM5ZNhZso9ZRiWWIhdE/sTmZsJnvL97LhwAY2FWzCZXcRExFDsacYp81JVlwWBZUFFFQWEO+Kx27ZsVk2fEEfgWCACp+cP6glgs0TqH8/bcqG2+HG4/e02N9kU7amAzKOd77PPgCl1KlAOfBiUwqgQdkfAb/WWp/e6onrKABbaRcCsTkt7+C3hyMOGrL2InH/nDqLgdFnsKr8s3qbpw+4nCuG/IzPcz6nU1wneiT24Kx/nUWiO5Gnz3uaKa9NkYJ7B0LaqtaqDgXdJVokrvWRKLUv6a6RkNVEmocmiKcrxRv7MK7zWXxe8JIIeAVR/ix8Xjt9usWzPn8j1UUpELO7+TBCbww2m+K87j/C5qokIyaDrNgsdpXuYlneMr7JlfwDLpsLb8Bb6x7oGNORrNgsthZt5UDlAdx2NylRKZzf83z6p/bnpVUvsXzfchyWg6AO4rQ5KasuI8IWgc2yU+wpauTyaExzLo7Q1hZ8v6HVQZu4BW2+1t0MuuaoSuG0RRCVdxaFH88AV7EI9fjtkLhNwu2i9kP68kahvy1ejXagg2BXEVgWRLukb8RhcxDjjKHCV0FZdRkef9utDpuyodHiGmg4/uEQsJSFVdN3EhKcqkYbHq6f3WE5cNjEpVFRXVHveA2fpV3ZiYmIwWbZCOogCa4ESrwllHnFWgoJf7fdTVxEHLZANNi9eINVWEEHnmAl3oCHsuqyRgpCoXBYDuw26atRShHUQXwBH267GwsH0c5InDYXygriD/rxB/04bU4cNof0LekgJeXVVO5Pxed1EJNSjMPlwWZX2CzweCwCldFU5nWRcQhZHjqnReFwBvFrH5VeLxXBInTATmVpJN7SGILeSDomxtKhUwVl3jJ8VgkFhVC2pyOuuGI0QRw6hoLSCrQnBpvDj7L5iXI7iE+pJCrOR4I7hupgNe9d8t736wJSSnUB3muDAvgPMF9r/Xyrx6yjAJyFg6hObHr+30PmyXUMvuLfrHA8SUZCMjk3b8LntTNihAx97zhoLUOfHUp1sKZFfqAHq65bw8DR+2H0YzDonxC9v+lj15jjja6pMhknsVS7dtX4i5vBGyORJw1/h+L6D4POMV3ZuUuh45tIWfl90lrERcAmeZn2DqFjlpfUbvsIEiArNotKj5/P161FlydBh9VN3hOJZLcIIALApt10tPclatsl9O6YQb61ksUf9STbcTLjr/yMFz//Cm/UNgLRORKq21JMvVYS0eNqWz+NTdmIckbhC/jwBX0ttlojHZForbGURZW/qtblkuROItGdSF55HtWBarx+b6ut24PpIA0J6JBQtFt2UqNSa/35XeO7io/aX0WELQLLsoiLiENrTZwrjkR3IvvK97E6fzWl3lJinDGU+8oprCpsg7I/NGxBF8HyVHT07vrjYqqjsAWi6GgN4bSup9A7uScje3TnlJ79WbncQffuMnJ23z6ZR2PBAhmFvHFj/QFf/ftLzqLQ4KyqKhkxHcqmGxUlg7e2bGlct8REGVhWURGeoMhmk091dfNZXUMkJMigtaYGeoUGkoHUJTQSPzQKubT0ex4I1hYFoJSKBHKB7s25f5RSM4AZAKQzLKQA3HmnU5U+r6ldDo86nYM3nXQT05L+zPDhkNEnl3c+3ctJfzup3p/1lE6n8OWfboT+L8NHj9LtjovY5pVk9T/u+2PeWNfETCUlWRDXwtwFDYghlbIPb4ONF4Avkv++kso/3l/HM7/vw9wNc/l6//t8tPUD8srzoDqSqEAnKtwbWjxmRkwGPRJ7sLFgI3nleaRFpVPsKcIT8GBhJxigxcFGsY4ksqrPYvveAioTvnNInOgAAB2YSURBVK0XBmspq/Ef3OcitfACRiWdQ1y3rbz8WiUBqxK6zpdBQKunwZbzZECL5YeMxaT0XYs//VuKdnaAwh5klE1m99aE2lwtY8dKqoXXXgtnn+zaFUp8BRS6F5GQtY/f39SV5M77KfWWMrXn5cRFOymsKgR/BPFRkViWYtMmGXZfN+nc6adL2oJVqyTXTlmZhqj9pI/4lnLnFsr2pMkztPxEqCgGnb6Fzn0OkByVQO/k3nirHAztOJgvP43hob+UoGNzuObaIIOHaGKtdDrQn7K8dFJTJZfLeef7sOJ2o7Vme9F2AgTIiMnA6/cyOG1ws2kjGnKg8gAr965k+d7lRNgiSHQnsrFgIwrFzFEziXfFs69iH/sr9rOvYh8FlQUMSR9CdkI2O0p2UOYtY2vRVoamD6VrfNcmzxtKdhYVJQLtYAnqIAcqD7CrZBc7Snawo3gHadFpdI7vjD/oJ3d3kNXL3cTbO5CdnkyJzuXFD9bz1cp9BK1KCb21e4l3x9M7qTcRDjvfrNpPtT0fZ1IeWX3z6J7chX1re1GwuTtj+mQTa09l+XeKtWvDOaUsS4R1U+khOnaU7w4dJFXDBx+EsmrK/l6vXL/TKfcjPj6cCmPjRjm2UrI+OVmOs2OH5PTxtpKhpu4o4ZaUQlMJF5vn2FQAPwUu01q3KXVWXQsgNvciSjPfbMtugjcKIhqY56Ud4cMn4Kc/bvp8KC7rMIt//XEIXDS9fvhn3hC6pSexjU9bPfWF6u88OnM8Vz/yFp8t2wajZhNX3Zfei+YTDMKSHWsATa++fjYuS+Wkm/5MUlQSEfknM/NnWZz72K+pTFgirgZtNelmUChO73o6K/NWc8CTT6e4Tiy/Zjm7SnaxYu8KnDYnSinyK/Ip9hTz9oa3KfWWUuYtJzkymU5xWVQHqnny3CfpFtub+x+ALzeu4etvqyAyX6wbZzmn9O3DL846jytOPas2p09yiqbf0BK+2LAKKpMZP6Qr83d+BFnfQOpqevimsvl1GcRWl3nzZMh+QYFk0ywqknz627fD66+HBfLTT8sf4eqrZdlub1uelpgY+ZP26yd/zCVLGieVS00NZ6fMyJCkafPnN328Z5+VPDgg+9x0k0zoEsLhkMRtw4ZJorlQ5saVKyXPy6oWPIV2u9Q1lHht+XJpgU6fLmm2o6KkBVkXraXF2pYEby3x2mvw6acimE45RXIHhXIPrVkj+X88HhFEJSWSmTN0LZdcIjOtuVzhbK//+58IzI4dJVXD5MmimH0+yUw5f75kAi0tlfNmZ4sC/sc/5BmFsm/WxbKkXN++IlCLimQGvP79w/c3JkZyBOXnSwt89WpJpZ2SIuUPJrdPQ+F6cMK2edxuuVc2myiFnj3FIti8We5dZaWcKzSPQ2Sk3KuqKnk+rv9v78yjq6quP/7dyQuTEARShBAgTJKKP6MQQURbFauAaKW2DogW1Dq1iANaqS4BW+rP2iX0h1hURLQLkaXiUChYZ5wQoxIcAsigCQYIAZVBEkje/v3xvZf73st7GV9yQ97+rPXWu/fcc+8999xz9z7j3q1YzrZs8fJk+3Z+Ky1bMr9LS3m99u1ZHtu1A559tmkqgBcAPKuqT9foxiEK4KjiM7G/c4wvtaasOx945mUs+M/nuPN3x2L7hdkcHP7mZ8DP7oueBggCwXYof/Ar3DihE+ZsmIzWw+Yhq3M/fLrdcSP54UQkvzMNFb3+A2S+jSHnrceWvevRt0M/vL/VsX07ez1XVDpccw3t9Xd1prmGua6UIND7VfS/4R6s38cxgS5tuyBJktCnQx9cN+g6HN/5eGR3yUbJjyVYkr8Ew7oPw4DOA2I++v79/ChGjGChW7Uqem0uGGThe+ABmtTdEDIhacIEHp89m4Xss89on8Vl1iya7+3UiYJk8WIKhRtvpFC//vrwe/34I4XRH/7A9LnzqN3/9HQ2zWfMoCBasiT8/Jo4VklNpeC59FIqmaps1iQl0StZcTEFf3Z25TglJXxXK1bw/VUl5CM580ya3B44kM57du3yTDy7ppFD+cUv2Br55hsKy8WLef9x42jR9JhjauaTGGD+HjpEC6G//334sZwcCvPHH69s1x/g++zTh8rnm5A1gpH2bUKpiwB1vbnVlkCA+de6NSsYRx/NNPfrR2Xm+lVYvJjfmAgrID170gRz5DOnprIVkJNDpX7xxaws/PgjlUpmJq1xbthABbVnD8vpvn1Uzn37Mk1dutAEdDRT2qHs2+f52HBbIvWl0R3CVKcARKQ9gC0AuqtqjUbOajUNtDwFCByCaDLn736bA3TLxbWdF+DR4vGME0wGZm1Blzbdsb2kFJjSDnj/dmDl3cB1A4G09VXeokXBSBz85DfAmglY+mYx5n7xFyydegMwaB4weE64+WWXXf2AuWtotyOEtm1pzXDFCgqnaPTrBzz5pGJ1bgUKvg7g9tv5YUUT3Hl5FCaXXFLZL2ykoHYZOZJCrrSUDmIizfseOkQztG+9RfvjrpetMWNYcxehGeVFi6jQ+vf3zn3/fZoXnjiR5xYV0bzx2rWs7fzpT17c7t354Uf67O3Z07MTD7Cmff/9rJmvXs2PrkcP1kDvvZeOUQYNoiXJ3FymAfCsSALhi3Xat+ezl5XRNvvTT7OGWVtWraKgDgap0M89l3m2eDGFdHY20/Pqqzx+6BDfx8MPexYj9+6l3+KVlReTxyQQoDnpwYO5v3Ila+vnnENh9e23VBiFheG17NRUmpdesyY8T1q0YJ7s3s24bdp4Qrw0Yly6ZUsKw4qKcGuarVrx/QaDvG6XLrRwum5duKDt1o0mpR95hBWO++/n+09PZyuoVStul5VR8fTrx/zcvJnl9IQTWJbS06lYU1PD07dpE/N7zhxWRgIBOjSaMoX3DaW0lOV67Vpa7uzRA0e817JGVQAisgjAGQDSAOwAMBVACgCo6lwnzngAI1Q1hriLct2aKoCCU4CXnwCCyRh9YSmWvl0EbDqXA3SlR9PAU7ttwITTgfyLgM/Gcp73+LOAxUuA/DH0v3vN0EqX7vvDNdiYdwzwsxleYKyZQHnjkN33GOT9exhaH8zEgYNldCJTGt1T96hRwLJldIjy738z7Oyz2Uzu1i22KdusLH4wF1zA/aefpnMZgMLg2GMprD/5xPPMFEnkqsWOHdlUj1QU773ndQG88gqFdWEha54zZrD2npxMQfbZZ1ROCxeGu6WcNo0f4s4Y4+Uuy5ZRgO/aBdxwQ2VhOGAAlQjAD/6uuygcJk/2WlIu33/P2hfA7omyMj5vSgoFzLBh7PooLOR1XnqJQu13v6Op3qGVi0K9yMtjevLyYndNtG9P5ZqSQv/F9bE+GYvIWntKCvcja+vR+qIHDeJ/QQEVS8eObKVs20aFM3o0W2l5EXM1kpOZn2edRWV48sns8nK7gzp1Yi34qaco4IcOpY+FHj34HlNTqWSKithidGvho0d73WFlZUzXc8+xHC5yLL+0aUMPcpdfHnvdQ3Ok2buEbH2gLw60dobepwVRIwMel14IZIW7s3rh9EJcdl4Gazitd2PaPzZjc4uX8HrJkzgw5z2U7ujOPuQW+4CMDzBgzDJ80c4xrby7D1CUA3TcCCx+HvihR6V0hH5IbdrQI1BRkees/JZbWGiLi+vW/J0+nbXA2tK5M4V4RgaFt9uVceWVwPjx7KoYN877kGLRvTsFe2QNEajaR+1RR7FbAmAN+JJL6AmrfXvWCFeu5AeemsoaXkEBa6yHDjFs/35PkPXuTe9Pquxq6t6dNfBPP/Xul5wMvPMOhU8gymLjDRtYc3a7OAIB9mXfdlvlc9xuqtDtYJBK7ssv2YpwvYV9+im9qrmmiX/6U9aGQz8pt4XSqxeFW2EhBWV2NgXj5s1UfOPHU9AmJ1MAq3JcJRj0ylho10tqKoXeCScwr1zzxuPGsTb8wQeseGzaxMrEeedRmCcnc6xj504q106d2JIZOdJ77tJSxkuJWEBeUUFPbAUFdFZTXMyKw/PPV+/cpUcPvutlyyp7hYtGSgrLabdu7JIJvf6VV7L8jBkT/X03d5q9Aui7eik2Dh5Nq5wf/aFG1+t9yhfYOuQyHOxA7wqzzp2Fm4ZMQpcuLKglJSzsAPDAA4qZM6WS05LLLwdunbEBg8Y/Day8C/8zIAW9e3M175lnhse9/npg7tzaPHF0brqJtZjWrelXN9L/LhDezeHSqhWFgSuE77iDAuTee1n7uvji6u/dogWVxdatFM779vEj79GDH3kskpN578iiE6kUBg5kTdf14JSW5n38/fuziywzk/uu440DB5iOOXPYdTN5cnTft8OGsUW0YAG7W2J1tbmUlDBf5s1jt87+kI7K9HQK6A0bmN7MTNZ+lyxhK6Q6h+z9+1OZZmdzEO+//6Xgzc9n3/+2bZ4Tlpde4lhNKD/+yApEJK4C+uIL5nfv3sCLL7J2PnJk5biVvMshXKE1FKps2b3xBvP4uONY0y8qojJJT+f7CgRY63/sMeZpmzYc4C0sZCvv6qt5rffeY0tjxQoqyBNPZOv5ggvYl+6XL+imQjwVAO1L+PBDVyimVfPLekFZJKr/bdyoev75qkjL11NGblJV1cce847Pnq2qqlpeHvsaLVuqvvaa6rBhqjNnqn77rWowqLp7d/i16vu77TbVwkLV1auZpooK/h86pHrrraopKVWfP3y4HmbPHo3KXXepZmWptm7tnZeczH+R8Osde6zq3r2qX32l2qlT5ft16KD6zjuqU6aEh998s+ry5arTpoWHp6WprljBdASDzPs2bbzj3bur5uZGT3cke/eq/vnPqldcwfewdKnq9OnMK/f6dWHPHtVrrmF6+vXj/09+opqeXvn5TzxRdf581fx81XffVf3LX1QnTlR98EHVoiKmsTrKylS//75uaU1kSkvr/o6bKwByNV5yOF4XqvWNa6IAMj44/BEGAlULxTfeUP3jH7k9YYLqiy+GH8/KUt28WfXLL6OfP2pU+P7DD1MBtG3L/RYtKp/TubPqTTdx+8MPVQsKVJ9/XjUzMzze2Wfz/7rrVEeO5Pb06XyZBw+qLlzIgj5mTOV7/Pzn/J83zwu7+mrVrVurLiRuXlT3CwRUt2zxzisqUn3iCSqZm29WnTrVi3vUUfwfOFB15crw+5WUqN54o+rbb0dPTzCoumQJr98UcRWKKtP4r3+pvvqqf+kxjFgkjgJo/7WmptZMkD3yCGudgOp997EGHC3ewoUaVhsGVLOzVTdsCI83YoTq3/9e+fzZs1njmzmTNcKTT1Y97rjoL+r888PPXb2aSgVQTUpS/etfVQcPjp7O7dtVu3XjdkYGWwnr16v27evFefZZ1dtv9wT4kiVUVFdc4QnrJ59U3blT9e67VW+5hTXYoUN57NJLVXfsqLqwVVSE50PHjqq7dlV9jmEYDUfiKIDkUj311JopgMmTmTnvvsva9FVXsTkfDMbuUpk6lQI0P5/nPvRQ9HgTJ3rboTXcQYMYdvfdsV/WVVcxzvDhXlN28eLYzzF8OFsfqmxNpKezy8Fl9erK57RsSSUWGf7RR7HTVVRUu6b12rWqd9zBrgzDMPwjngqgaQ0C/9DNM7D23u3Aq3/D6adzhkcoM2Zwel8ov/wlB8hcRo3i4G9uLjBkCOeWR/Laa5y5E8rChRzQmjSJ+5mZHMR84glOOXzqKQ7Ylpd7syTy8qLPxXepqAgfuKqo4MrHmTN5rbw8Oqu+7TYOllVHaSlncrz1Fgf9/vY3hot4M3sCAeCii6q/lmEYRxbNdxbQzC3ALc6Kpf/lPP/jj+dij1C++85bEDV/PueV9+njzSMHuGAoI4NT4SZP5jTEUK67jqteI6e6AZy1kZ7Oe6xb5y3OCgY5u2bcOE7DGz48+irYmuJO6XOX6teVTz7hDJHU1PpfyzCMpk08FUCTEBdJh9px44dML7CMYZHC/6KLOGXxxhuB++6jCYNJkzidLHQRzrZt3gKiCy+sfM+5cyks10Wxs9a1K1sd33wTvjI3P59znseM8VoONZluGYukpPgI7IEDuUTehL9hGLWhSSyjCO7I8ixqzn8HOHYpIo2NuTznGOScM8cLO+kkztPOy+OKxvJydv+4CuC007y4CxZw8Y9rDRCgNcSePcPvE3oOwC6kW28ND8vMrGzUyzAM40ihSSgAtNvm2dMpOI2/EGKZPHBxl/Z/9BEVwI4dHAYNNSEwfjyFf1ISBXeoMHeNSsVi2TLg/PN5zdGjuQCpooKmBQzDMI5UmogCKAJ29Y95+IQTqlYAGRnsy//6a+67K3xDFcDjjwPLlwN3383VmQ895B0bO5bKQ5UDs8Egu3969aKVybFjqYTuuYfGwOpiN90wDKOp0TR6jZOCQHlsu6rRlsmHkpzMGvzq1eyjj6YAkpI48FtQwIFbgGMDrVtze+ZMmjxeu5b2XXr3ptK4+GJ2Kf3mNzQ0ZsLfMIzmQtNoAQBARbgCCLUfnpLCQddoNtxdMjM5rfPNN73++0grkhMn0lZOSQltxyxaRDs4/ft7jkVC7+E6LQGAO++s22MZhmE0VZqOAohoAbRqxdk2nTrRfOwxx1R9eqi9+3ffZasg0gFDWprngWrIEIZlZNAwWDAY20JmYWFlm+SGYRhHOk1HAUS0AA4epNAN7auvih49wvdzcqJbDezQgaZwXUcbLklJtL3/wQccS3jmGSqfPXuoJAzDMJobTWMMAKjUAigvr7lLPIAmm++5x9tfsCB6vOXLaZ8+J8YyinHj+J+TQxO0v/pVzdNgGIZxJOG/Agg6SQiEex1x3dbVlLZt6UDl/fe5MjYrK3q8QYNYu3fXAEQyfToXh/XpU/N7G4ZhHIn43wWU5Lo48nwkrl3LqZ+1aQG41NfdX1JSuP9bwzCM5kq1LQARmS8ixSLyeRVxzhCRNSLyhYi8XaeUrJp0eNM16VAXBWAYhmHUjJq0ABYAeAhAFGeFgIgcDeBh0Cl8gYjUfqb8vQeBoGeVzRSAYRhGw1NtC0BVVwLYXUWUsQCWqGqBE7+41qkIhpvkdBVAbcYADMMwjNoRj0HgYwF0EJG3RORjEbkyVkQRuVZEckUk93BgsHISTj2V/9YCMAzDaDjiMQgcADAIwHAArQF8ICKrVHVDZERVfRTAo4DjDwAA1oyPeWFTAIZhGA1HPBTAVgAlqrofwH4RWQkgG0AlBRCdYMwj1gVkGIbRcMSjC+glAKeLSEBE2gAYAiC/xmcnl8c8ZC0AwzCMhqPaFoCILAJwBoA0EdkKYCqAFABQ1bmqmi8iKwCsBavz81Q15pTRSpTEnnRvCsAwDKPhqFYBqOplNYjzAIAH6pSCdVH8NTpYF5BhGEbD0QRMQXhTQHv3Dj9kLQDDMIyGw38FEOL7d+zY8EOmAAzDMBoO/xVA0FMAoSade/SIbs7ZMAzDiA9NQAF4wxChTl9WrfIhLYZhGAmE/wogpAuofXsvuGVsF8GGYRhGHPBfAYS0ADp29IJj2es3DMMw4kMTUABsAQQCwIABXrC1AAzDMBoW/xWA0wWUlUVnLC4B/13VGIZhNGv8F7NOF5Db5bN8OfDCC4CIj2kyDMNIAJqAAmALwJ3yOWIEf4ZhGEbD4n8XkNMCsDn/hmEYjYv/CkC9QWDDMAyj8WgCCoBJsBaAYRhG4+K/AnCwFoBhGEbj0mQUgLUADMMwGhf/FECIGWjAWgCGYRiNjX8KYOdxwGMfHt61FoBhGEbj4mMLIAB8O/jwri38MgzDaFyqVQAiMl9EikUkqp9fETlDRH4QkTXO7566JKSsrC5nGYZhGHWlJj3vCwA8BOCpKuK8o6qj65IAEUAVOHCgLmcbhmEYdaXaFoCqrgSwu6ES4Hb9mAIwDMNoXOI1BjBURPJEZLmIDKg+uof1/RuGYfhDPBTAJwB6qmo2gNkAXowVUUSuFZFcEcmNPNa5cxxSYhiGYdSYeisAVd2jqvuc7f8ASBGRtBhxH1XVHFXNccPcFsDQofVNiWEYhlEb6r38SkS6ANihqioig0Glsqum57doASxdCpx9dn1TYhiGYdSGahWAiCwCcAaANBHZCmAqgBQAUNW5AH4N4AYRKQdwAMClqqo1TYAIcO65dUi5YRiGUS+kFrI6vjeWHAVy0bYtsHevL0kwDMM44hCRj0O70euD78bgbBaQYRiGP5gCMAzDSFBMARiGYSQopgAMwzASFFMAhmEYCYopAMMwjATFdwVgGIZh+IPvCsBaAIZhGP5gCsAwDCNBMQVgGIaRoPiuAFJS/E6BYRhGYuK7AujWze8UGIZhJCa+KwBrARiGYfiD7wogyfcUGIZhJCa+i9/kZL9TYBiGkZiYAjAMw0hQTAEYhmEkKL4rABsDMAzD8Affxe8rr/idAsMwjMSkWgUgIvNFpFhEPq8m3skiUiEiv45f8gzDMIyGoiYtgAUARlQVQUSSAdwPwOrzhmEYRwjVKgBVXQlgdzXRJgJ4HkBxPBJlGIZhNDyB+l5ARLoBGAPgLAAnVxP3WgDXOrtlALuVzCAc0gCU+J2IJoLlhYflhYflhUf/eF2o3goAwCwAf1TVCqlGkqvqowAeBQARyVXVnDjc/4jH8sLD8sLD8sLD8sJDRHLjda14KIAcAM84wj8NwCgRKVfVF+NwbcMwDKOBqLcCUNVe7raILACw1IS/YRhG06daBSAiiwCcASBNRLYCmAogBQBUdW497v1oPc5tblheeFheeFheeFheeMQtL0RV43UtwzAM4wjC95XAhmEYhj+YAjAMw0hQfFEAIjJCRNaLyEYRudOPNDQWItJdRN4UkXwR+UJEJjnhHUXkVRH5yvnv4ISLiPyfkzdrRWSgv08Qf0QkWUQ+FZGlzn4vEfnQyYvFItLCCW/p7G90jmf6me54IyJHi8hzIrLOKR9DE7VciMgtzvfxuYgsEpFWiVQuopncqUtZEJHfOvG/EpHfVnffRlcAjtmIOQBGAjgOwGUiclxjp6MRKQdwm6r+FMApAH7vPO+dAF5X1X4AXnf2AeZLP+d3LYB/Nn6SG5xJAPJD9u8HMNPJi+8AXO2EXw3gO1XtC2CmE6858Q8AK1Q1C0A2mCcJVy6cxaQ3AchR1eMBJAO4FIlVLhagssmdWpUFEekITtIZAmAwgKmu0oiJqjbqD8BQAK+E7E8BMKWx0+HXD8BLAH4BYD2Ark5YVwDrne1HAFwWEv9wvObwA5DhFOazACwFIOAKz0Bk+QBtSw11tgNOPPH7GeKUD6kAtkQ+TyKWCwDdABQC6Oi856UAzk20cgEgE8DndS0LAC4D8EhIeFi8aD8/uoDcl+2y1Qlr9jhN1ZMAfAjgGFXdBgDOf2cnWnPPn1kA7gAQdPY7AfheVcud/dDnPZwXzvEfnPjNgd4AdgJ4wukOmyciRyEBy4Wqfgvg7wAKAGwD3/PHSMxyEUpty0Kty4gfCiCavYhmPxdVRNqCBvNuVtU9VUWNEtYs8kdERgMoVtWPQ4OjRNUaHDvSCQAYCOCfqnoSgP3wmvjRaLZ54XRT/BJALwDpAI4CuzkiSYRyURNiPX+t88UPBbAVQPeQ/QwART6ko9EQkRRQ+C9U1SVO8A4R6eoc7wrPkmpzzp9hAC4Qka8BPAN2A80CcLSIuIsSQ5/3cF44x9ujesu0RwpbAWxV1Q+d/edAhZCI5eJsAFtUdaeqHgKwBMCpSMxyEUpty0Kty4gfCuAjAP2cEf4W4GDPyz6ko1EQGkl6HEC+qj4YcuhlAO4o/W/BsQE3/EpnpP8UAD+4zcAjHVWdoqoZqpoJvvc3VPVyAG8CcB0JReaFm0e/duI3i5qeqm4HUCgirmXH4QC+RAKWC7Dr5xQRaeN8L25eJFy5iKC2ZeEVAOeISAenVXUOqvPR4tNgxygAGwBsAnCX34MvDfysp4HNsLUA1ji/UWCf5esAvnL+OzrxBZwltQnAZ+DMCN+fowHy5QzQbhTA/vDVADYCeBZASye8lbO/0Tne2+90xzkPTgSQ65SNFwF0SNRyAWA6gHUAPgfwLwAtE6lcAFgEjn8cAmvyV9elLAC4ysmXjQAmVHdfMwVhGIaRoNhKYMMwjATFFIBhGEaCYgrAMAwjQTEFYBiGkaCYAjAMw0hQTAEYhmEkKKYADMMwEpT/BwWPsCPI/5zwAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "########################\n", "## Plot 10 trajectories\n", "########################\n", "\n", "##############################################################\n", "# Simulate M=10 trajectories of the estimators I_n et I_n^c\n", "# Wished size for the arrays: M x N\n", "##############################################################\n", "N = 1000\n", "integers1toN = np.arange(1, N+1)\n", "\n", "M = 10\n", "\n", "X = np.random.uniform(low=0, high=1, size=(M,N))\n", "\n", "Y = np.exp(X)\n", "Y_c = 1 + X\n", "\n", "I = np.cumsum(Y, axis=1) / integers1toN\n", "I_c = m_c + np.cumsum(Y - Y_c, axis=1) / integers1toN\n", "\n", "############\n", "## Plotting\n", "############\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "\n", "ax.plot(integers1toN, I[0], color=\"b\", label=\"MC\")\n", "ax.plot(integers1toN, I[1:].T, color=\"b\")\n", "\n", "ax.plot(integers1toN, I_c[0], color=\"g\", label=\"Control\")\n", "ax.plot(integers1toN, I_c[1:].T, color=\"g\")\n", "ax.axhline(Exp_gY, color=\"r\", label=\"True expectation\")\n", "\n", "ax.set_xlim(0, N)\n", "ax.set_ylim(1.4, 2.2)\n", "ax.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 2 (b):\n", "\n", "Plot the histograms of the (unnormalized) errors\n", "\n", "$$\\bigl(I^j_N - \\mathbb{E}[Y] \\bigr)_{1 \\le j \\le M} \\quad \\mbox{and}\n", "\\qquad \\bigl(I^{c,j}_N - \\mathbb{E}[Y] \\bigr)_{1 \\le j \\le M}$$\n", "\n", "where $\\bigl(I^j_N\\bigr)_{1\\leq j\\leq M}$ and $\\bigl(I^{c,j}_N\\bigr)_{1\\leq j\\leq M}$ are i.i.d. samples of the two estimators $I$ and $I^c$, for a fixed value of $N$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 500 \n", "M = 1000 \n", "\n", "####################################################\n", "# The exact value of the expectation, for comparison\n", "Esp_gY = np.exp(1.)-1.\n", "\n", "############################################\n", "# The exact value of m_c, again\n", "m_c = 1.5\n", "\n", "############################################\n", "# Simulate the required samples of Y and Y_c\n", "Y = ????\n", "Y_c = ????\n", "############################################\n", "\n", "############################################\n", "# Construct the samples of size M for the\n", "# two estimators, without and with control variate\n", "I_N = ?????\n", "Ic_N = ????\n", "\n", "#########################################\n", "## Plot the histograms of the errors\n", "\n", "plt.hist( ????? , density=\"True\", bins=int(np.sqrt(M)), label=\"MC\")\n", "plt.hist( ????? , density=\"True\", bins=int(np.sqrt(M)), label=\"Control\")\n", "plt.title(\"Estimation with and without control variate, N = %1.0f\" %N)\n", "\n", "plt.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 3: Optimal control variate\n", "\n", "We now introduce a parameter $\\lambda$ and consider the estimator \n", "\n", "\\begin{eqnarray*}\n", "I_n^\\lambda = \\lambda \\, m_r + \\frac{1}{n}\\sum_{i=1}^n \\bigl(f(X_i) - \\lambda \\, f_c(X_i)\\bigr),\n", "\\qquad \\lambda \\in \\mathbb R.\n", "\\end{eqnarray*}\n", "\n", "Propose a choice for the parameter $\\lambda$.\n", "\n", "Plot the trajectoires of the resulting estimator $I_n^\\lambda$ and compare them with $I_n$, then compare the histograms of the errors for a fixed value of $n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "##############################################\n", "## Empirical estimation of the optimal lambda\n", "##############################################\n", "n_1 = 100\n", "X = np.random.rand(n_1)\n", "Y = np.exp(X)\n", "\n", "## Empirical optimal lambda\n", "lambda_opt = ???\n", "\n", "print(\"Estimation of optimal lambda = %1.3f\" %lambda_opt)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#######################################\n", "### For the histograms: M x N samples\n", "#######################################\n", "M = 1000\n", "N = 500\n", "integers1toN = np.arange(1,N+1)\n", "\n", "Y = ?????\n", "Y_control = ?????\n", "\n", "#####################################\n", "### Plot the first 10 trajectories\n", "#####################################\n", "I_n = np.cumsum(Y[0:10,:], axis=1) / integers1toN\n", "Ic_n = np.cumsum(Y_control[0:10,:], axis=1) / integers1toN\n", "\n", "## Plot the trajectories\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "\n", "ax.plot(integers1toN, I_n[0], color=\"b\", label=\"MC\")\n", "ax.plot(integers1toN, I_n[1:].T, color=\"b\")\n", "\n", "ax.plot(integers1toN, Ic_n[0], color=\"g\", label=\"Control\")\n", "ax.plot(integers1toN, Ic_n[1:].T, color=\"g\")\n", "ax.axhline(Esp_gY, color=\"r\", label=\"True expectation\")\n", "\n", "ax.set_xlim(0, N)\n", "ax.set_ylim(1.4, 2.2)\n", "ax.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#############################################\n", "## For the histograms: evaluate the errors\n", "## of the two estimators I_N et I_N control\n", "#############################################\n", "error_N = ?????\n", "error_control_N = ?????\n", "\n", "plt.hist( error_N , density=\"True\", bins=int(np.sqrt(M)), label=\"MC\")\n", "\n", "plt.hist( error_control_N, density=\"True\", bins=int(np.sqrt(M)), label=\"Controle\")\n", "\n", "plt.title(\"Estimation with and without control variate, N = %1.0f\" %N)\n", "\n", "plt.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2. Importance sampling\n", "\n", "We consider input random variables $Y$ with standard Gaussian distribution $\\mathcal N(0,1)$.\n", "\n", "The aim is to evaluate $\\mathbb E [g(Y)] = \\mathbb E \\bigl[ (Y - 2)^+ \\bigr]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 1: standard Monte-Carlo\n", "\n", "Simulate and plot the trajectories of the empirical mean estimator $I_n = \\frac 1 n \\sum_{i=1}^n (Y_i - 2)^+$ obtained from $n$ i.i.d copies $(Y_i)_{1 \\le i \\le n}$ of $Y$, then plot the histogram of the errors." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Function g: the positive part (x - 2)^+\n", "def g(x):\n", " return np.maximum(x-2.,0.)\n", "\n", "N = 2000\n", "integers1toN = np.arange(1,N+1) \n", "\n", "############################################\n", "# Complete with N draws of the standard\n", "# Gaussian distribution\n", "Y = ?????\n", "\n", "# Evaluate the function g(Y) on the sample\n", "GY = ?????\n", "\n", "#############################################################\n", "# Stock into the variable 'mean' the MC estimation E[g(Y)],\n", "# into 'var' the empirical variance, and into 'radius_CI'\n", "# the half-lenght of the 95% confidence interval\n", "mean = ?????\n", "var = ?????\n", "radius_conf_int = ?????\n", "\n", "print(\"MC estimation = %1.4f\" %mean )\n", "print(\"95%% confidence interval for E[g(Y)] = [ %1.4f , %1.4f ] \\n\" \\\n", " %(mean - radius_conf_int, mean + radius_conf_int))\n", "print(\"Relative error = %1.2f\" %(radius_conf_int / mean))\n", "\n", "######################################\n", "# Trajectories of the empirical mean\n", "######################################\n", "M = 10 # Number of trajectoires\n", "\n", "#####################################################\n", "# Evaluate the M trajectories of the empirical mean\n", "# I_n for n ranging from 1 to N\n", "I_n = ?????\n", "\n", "##########################\n", "# Plot the trajectories\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "ax.plot(integers1toN, I_n.T, color=\"b\")\n", "\n", "ax.set_xlim(0, N)\n", "ax.set_ylim(0, 3.e-2)\n", "ax.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 2: Importance sampling\n", "\n", "**(a) (theoretical)** Show that, if $\\tilde Y$ is a Gaussian random variable centered at $\\theta \\in \\mathbb R$ and with unit variance, then\n", "\n", "$$\n", "\\mathbb{E} [ g(Y) ]\n", "=\n", "\\mathbb{E} \\Bigl[ g(\\tilde Y) \\, e^{-\\theta \\, \\tilde Y + \\frac{\\theta^2}2} \\Bigr],\n", "\\qquad \\theta \\in \\mathbb R\n", "$$\n", "\n", "\n", "What is the interest of such a formula?\n", "\n", "**(b)** \n", "Construct an estimator $J_n$ of $\\mathbb{E}[g(Y)]$ based on the simulation of the gaussian distribution with unit variance and centered at $\\theta=2$.\n", " \n", "Plot the trajectories of the estimator $J_n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "theta = 2.\n", "\n", "N = 2000 # Taille echantillon\n", "integers1toN = np.arange(1,N+1)\n", "\n", "###############################################################################\n", "# Complete with N draws of the Gaussian distribution centered at theta = 2\n", "# and evaluate the importance sampling estimator constructed with these N samples\n", "###############################################################################\n", "J_N = ????\n", "\n", "# Empirical variance and confidence interval\n", "var_IS = ????\n", "radius_conf_int_IS = ????\n", "\n", "############################################\n", "# The gain in terms of variance\n", "variance_ratio = var_IS / var\n", "############################################\n", "\n", "print(\"Importance sampl estimator = %1.4f\" %J_N)\n", "print(\"95%% confidence interval for E[g(Y)] = [ %1.6f , %1.6f ]\" \\ \n", " %(J_N - radius_conf_int_IS, J_N + radius_conf_int_IS))\n", "print(\"Relative error = %1.4f\" %(radius_conf_int_IS / J_N))\n", "\n", "####################################################\n", "# Trajectories of the importance sampling estimator\n", "####################################################\n", "M = 10\n", "\n", "#######################################################\n", "# Complete with the M trajectories of the IS estimator\n", "# J_n, for n ranging from 1 to N\n", "J_n = ?????\n", "\n", "# Plot the M trajectories\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "ax.plot(integers1toN, J_n.T, color=\"b\")\n", "\n", "ax.set_xlim(0, N)\n", "ax.set_ylim(0, 3.e-2)\n", "ax.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 3: \n", "Compare the histograms of the errors of the estimators $I_n$ et $J_n$, for a fixed value of $n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 1000 \n", "M = 1000 # Sample size for the histogram of the error\n", "\n", "############################################\n", "# Complete with M x N draws \n", "# + from the distribution N(0,1)\n", "# + and from the distribution N(2,1)\n", "# and evalute the estimators\n", "GY = ?????\n", "\n", "GY_importance = ?????\n", "\n", "############################################\n", "# Evaluate the errors for the two estimators\n", "# Wished output size: sample of size M\n", "error_MC = ?????\n", "error_Importance = ?????\n", "\n", "# Plot the histograms of the errors\n", "plt.hist( ???, normed=\"True\", bins=int(np.sqrt(M)), label=\"error MC\")\n", "\n", "plt.hist( ???, normed=\"True\", bins=int(np.sqrt(M)), label=\"error Imp sampling\")\n", "\n", "plt.legend(loc=\"best\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ In practice:\n", "\n", "The optimal value $\\theta^*$ is the one that minimizes the variance of the importance sampling estimator: \n", "\n", "$$\n", "V(\\theta) = \n", "\\mathbb E \\left[ g(Y + \\theta) \\, e^{-2 \\, \\theta \\, Y - \\, \\theta^2} \\right]\n", "- \\mathbb E [g(Y)]^2.\n", "$$\n", "\n", "In practice, one can look for a value of $\\theta$ that minimizes an empirical approximation $V_n(\\theta)$ of $V(\\theta)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Application: estimation of VaR" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\newcommand{\\R}{\\mathbb{R}}$\n", "$\\newcommand{\\VR}{\\operatorname{VaR}}$\n", "Let $(X_i)_{i\\ge 1}$ be a sequence of i.i.d. random variables with density $f$, assumed to be continuous and strictly positive on an interval $I$ (and zero outside $I$, if $I$ is bounded).\n", "The related _cumulative distribution function_ (cdf) is\n", "\n", "$$\n", "F(x) = \\mathbb{P}(X_1\\leq x) = \\int_{-\\infty}^x f(u) \\textrm{d} u, \\qquad x \\in \\R.\n", "$$\n", "\n", "The function $F:\\R \\rightarrow [0,1]$ being continuous and strictly increasing on $I$ (under the hypotheses above on the density $f$), the quantile function $Q : [0,1] \\rightarrow \\R \\cup \\{\\pm \\infty\\} $ is defined by:\n", "\n", "$$\n", "Q(u)=F^{-1}(u) \\quad \\textrm{ for all } \\quad u \\in (0,1), \n", "\\qquad Q(0)=\\inf I,\n", "\\qquad Q(1)=\\sup I.\n", "$$\n", "\n", "The goal of this exercise is to estimate the quantile function $Q$ from an i.i.d. sample $X_1, \\ldots, X_n$.\n", "When $X$ represents the loss of a position, the $\\operatorname{VaR}$ is defined by\n", "\n", "$$\n", "\\VR(\\alpha)=Q(\\alpha) \\qquad\n", "\\quad 0<\\alpha<1,\n", "$$\n", "\n", "for values of $\\alpha$ close to $1$.\n", "\n", "We will consider $X_i$ with a standard Gaussian distribution (the functions `norm.cdf` and `norm.ppf` from `scipy.stats` can be used to compute the exact values of $F$ et $Q$).\n", "\n", "For every $n$, we denote $(X_{(i)})_{1\\leq i\\leq n}$ the variables $(X_i)_{1\\leq i\\leq n}$ sorted in the increasing order: \n", " \n", "$$\n", "X_{(1)} Q_n(u)\n", "########################################################\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\newcommand{\\maj}{>}$\n", "#### Question 3: estimation of quantiles by Importance sampling \n", "\n", "The application of the importance sampling (IS) technique to the estimation of quantiles (hence of VaR) can be implemented in different ways. \n", "The approach we follow below amounts to replacing the standard empirical approximation $F_n$ of the cdf of $X$ with the approximation obtained from IS.\n", "\n", "Since we are interested in the right tail of the distribution of $X$ (representing a loss value), we rather focus on the survival function $\\overline F(x) = \\mathbb P(X > x)$.\n", "It is immediate that\n", "\n", "$$\n", "Q(u) = \\inf\\{x \\in \\R : F(x) \\ge u\\}\n", "= \\inf\\{x \\in \\R : \\overline F(x) \\le 1- u\\},\n", "\\qquad u \\in (0,1).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The central point of the method is the following: instead of the standard Monte-Carlo approximation $\\overline F_n(x) = \\frac 1n \\sum_{i=1}^n 1_{X_i \\maj x}$ of the survival function, we consider the approximation given by the IS estimator:\n", "\n", "$$\n", "\\overline F_n^{\\theta}(x) =\n", "\\frac 1n \\sum_{i=1}^n\n", "1_{X_i + \\theta \\, \\maj x}\n", "\\,\n", "L(X_i),\n", "\\qquad \\mbox{where }\n", "L(G_i) = \\exp \\Bigl( -\\theta \\, X_i - \\frac{\\theta^2}2 \\Bigr).\n", "$$\n", "\n", "Note that, exactly as $\\overline F_n$, the function $\\overline F_n^{\\theta}$ is also piece-wise constant between the points $X_i$.\n", "The generalized inverse of $\\overline F_n^\\theta$, then, will also take its values inside the set of points $(X_i)_{1 \\le i \\le n}$.\n", "\n", "\n", "If we now define the _quantile by Importance sampling_ $Q_n^{\\theta}(u)$ by\n", "\n", "$$\n", "Q_n^{\\theta}(u)\n", "=\n", "\\inf\\{\n", "x \\in \\R:\n", "\\overline F_n^{\\theta}(x) \\le 1-u\n", "\\},\n", "$$\n", "\n", "we have\n", "\n", "$$\n", "Q_n^{\\theta}(u) = X_{(i(u))} + \\theta,\n", "\\qquad\n", "\\mbox{where } i(u) =\n", "\\max \\Bigl\\{ k : \\frac 1n \\sum_{j = k}^n L(X_{(j)})\n", "\\le (1 - u) \\Bigr\\}\n", "$$\n", "\n", "which is the estimator we want to implement." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Quantile level\n", "u = 0.95\n", "Var_u = sps.norm.ppf(u)\n", "\n", "# Simulations\n", "n = int(1.e3)\n", "\n", "# iid sample of the Gaussian distribution N(0,1), size=n\n", "X = np.random.randn(n)\n", "\n", "#############################################\n", "# The indexes that allow to sort the array X\n", "# in the increasing order \n", "indexes = np.argsort(X)\n", "\n", "################################\n", "## Standard empirical quantile\n", "################################\n", "sorted_X = X[indexes]\n", "\n", "empirical_VaR = sorted_X[ int(n*u) ]\n", "\n", "####################\n", "## *Quantile by IS*\n", "####################\n", "\n", "############################################\n", "# A possible choice for the parameter theta \n", "theta = empirical_VaR\n", "\n", "weights_L = np.exp(-theta * sorted_X - (theta**2)/2)\n", "\n", "###############################################################\n", "# In order to compute the cumulated sums backwards from n to k,\n", "# we reverse the array of weights L(G)\n", "weights_L_inverted = weights_L[::-1]\n", "\n", "cumulated_sum_weights = np.cumsum(weights_L_inverted)\n", "\n", "##################################################################\n", "# Complete with the evalution of the empirical quantile Q^{theta}_n \n", "# with IS, as defined above\n", "##################################################################\n", "empirical_VaR_IS = ???\n", " \n", "print(\"u = %1.2e, n = %1.0e \\n\" %(u, n) )\n", "\n", "print(\"True value of the VaR = %1.3f \\n\" %Var_u)\n", "\n", "print(\"Empirical VaR = %1.3f \\n\" %empirical_VaR)\n", "\n", "print(\"Empirical VaR by IS = %1.3f \\n\" %empirical_VaR_IS)\n", "\n", "\n", "################################\n", "# We can also plot the estimated\n", "# survival function: \n", "# Remove the if 0\n", "################################\n", "\n", "if 0:\n", " ###############################\n", " ## We fix some values for the x axis\n", " x_grid = sorted_X[n-80 : ]\n", "\n", " theoretical_values = sps.norm.sf(x_grid)\n", "\n", " plt.plot(x_grid, theoretical_values, \"k+\", label=\"True survival function\", markersize=13)\n", "\n", " ###############################\n", " # Standard MC estimate \n", " MC_estimates = np.arange(n, 0, -1, dtype=float) / n\n", "\n", " plt.plot(x_grid, MC_estimates[n-80:], color=\"red\", label=\"MC estimation\", linewidth=1.5)\n", "\n", " ###############################\n", " # Estimation by IS\n", " IS_estimates = np.zeros(x_grid.size)\n", "\n", " for i, x in enumerate(x_grid):\n", " IS_estimates[i] = np.mean( (X + theta > x) * np.exp(-theta * X -(theta**2)/2) )\n", "\n", " plt.plot(x_grid, IS_estimates, color=\"blue\", label=\"Estimation by IS\", linewidth=1.5)\n", "\n", " plt.title(r\"Estimation of the survival function $P(X > x)$\", fontsize=15)\n", " plt.xlabel(\"x\", fontsize=12)\n", " plt.legend(loc=\"best\", fontsize=12)" ] } ], "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }