py_05_simulations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 4:  Simulation and Bootstrap\n",
    "## ECON5170 Computational Methods in Economics\n",
    "#### Author: Zhentao Shi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\zhent\\AppData\\Local\\Temp\\ipykernel_20464\\1965519564.py:2: DeprecationWarning: \n",
      "Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),\n",
      "(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)\n",
      "but was not found to be installed on your system.\n",
      "If this would cause problems for you,\n",
      "please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466\n",
      "        \n",
      "  import pandas as pd\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "random.seed(888)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simulation\n",
    "\n",
    "\n",
    "* check finite-sample performance\n",
    "* bootstrap, a data-driven inference procedure\n",
    "* generate non-standard distributions\n",
    "* approximate integrals with no analytic expressions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Appetizer\n",
    "\n",
    "calculating $\\pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAHFCAYAAAAwicqeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB7zElEQVR4nO3de1xM+f8H8NfMVBOlC1FCLotiXSJqw8olhBVrsSWRclmXXda1dheLtbntza7FWsqd3HO/JOWSkLts5R4qhKlQqvn8/phf8zW618ycc2bez8djHqOZM+e856h5z/lc3h8RY4yBEEIIIbwm5joAQgghhJSOEjYhhBAiAJSwCSGEEAGghE0IIYQIACVsQgghRAAoYRNCCCECQAmbEEIIEQBK2IQQQogAUMImhBBCBIASNiGEECIAgkrY0dHR6NevH2xtbSESibBnz55SX3Py5Em0bdsWUqkUjRs3RmhoaKFtli9fjgYNGsDY2BguLi44f/68+oMnhBBCKkFQCfv169do3bo1li9fXqbt7927h759+6Jr1664cuUKJk+ejFGjRuHIkSPKbbZt24YpU6Zgzpw5uHTpElq3bo1evXrh6dOnmnobhBBCSLmJhLr4h0gkwu7duzFgwIBit5k5cyYOHDiAGzduKB/z8vLCq1evcPjwYQCAi4sL2rdvj7/++gsAIJfLUa9ePXz99dcIDAzU6HsghBBCysqA6wA0KSYmBu7u7iqP9erVC5MnTwYAvHv3DnFxcQgKClI+LxaL4e7ujpiYmGL3m5OTg5ycHOXPcrkcL168QI0aNSASidT7JgghhPAaYwyZmZmwtbWFWKy5hmudTtipqamwtrZWecza2hoZGRl4+/YtXr58ifz8/CK3+e+//4rdb3BwMObOnauRmAkhhAhTcnIy6tatq7H963TC1pSgoCBMmTJF+bNMJoOdnR3++OMPODo6chcY0bp58+xw4oQlGCvcsiISMXTr9hKzZz/U2PGzs7Px4MEDPHz4EMnJyXj06BEePnyIR48e4e3bt/8fhwg2NjaoWbMmatSooXKzsrJS3puYmHDSQpSTk4P09HQ8f/4c6enphf6dnp6OJ0+e4N27dwAUrWB16tRBvXr1ULduXdjZ2aFu3bqoX78+qlevrvX4K+LQIUssWmQHkQhgDMr7mTMfonfvl1yHR8rpypUrmDRpEqpVq6bR4+h0wraxsUFaWprKY2lpaTAzM0OVKlUgkUggkUiK3MbGxqbY/UqlUkil0kKPOzo6onPnzuoJnghC+/bAyZNAfn7h58RioH17c3Tu3EAtx8rKysKVK1cQFxenvN26dQsFw1Csra1hb28PNzc3NG3aFPb29mjatCkaNWoEIyMjtcTAFblcjsePHyMxMVF5S0hIQFxcHHbs2AG5XA4AsLW1hZOTE5ycnNCuXTs4OTmV+LfMhaQkYPFiRYIuGEFUcL94cUsEBACNG3MXH6k4TX/h1emE7erqioMHD6o8duzYMbi6ugIAjIyM4OTkhIiICOXgNblcjoiICEycOFHb4RIB8vdXfPgWhTEgIKBi+83NzcXFixdx/vz5QslZKpWiVatW6Ny5M6ZMmYJWrVqhadOmMDc3r/gb4TmxWIx69eqhXr166N69u8pz7969w507d3Dz5k3luVq2bBlevHgBoHAS79ixIywsLDh4Fwpr1yquqIsiEgFr1gDBwdqNiQiDoBJ2VlYWbt++rfz53r17uHLlCqpXrw47OzsEBQXh8ePHWL9+PQDgq6++wl9//YUZM2bA398fJ06cQFhYGA4cOKDcx5QpUzBixAi0a9cOzs7O+P333/H69WuMHDlS6++PCE+TJooP2IAAFGreXLOm7FdKubm5iIuLQ2RkJE6ePInTp0/jzZs3kEqlaN26Ndzc3DBlyhQ4OTnh448/hqGhoWbfmIAYGRmhWbNmaNasGQYNGgRAMQjowYMHuHjxYqEkLhKJ0KZNG3Tt2hVdunTBp59+qtUvO/fv/++K+kOMKZ4npEhMQCIjIxmAQrcRI0YwxhgbMWIEc3NzK/QaR0dHZmRkxBo1asRCQkIK7ffPP/9kdnZ2zMjIiDk7O7Nz586VKy6ZTMYAsKioqAq+MyJ0SUmMBQYy5uWluE9KKnn73Nxcdu7cObZw4ULWq1cvZmJiwgAwU1NT1rt3b7Zo0SIWGxvL3r17p503oAfkcjm7ffs2W7NmDRs2bBirU6cOA8DEYjFr164dmz59Ojtw4ACTyWQajSMwkDGJpKBBXPUmkSieJ8ISFRXFAGj8d0ew87D5JCMjA+bm5oiKiqI+bFKsFy9e4NChQ9i3bx8OHTqEjIwMmJiYoFOnTsqrPScnJxgYCKrhS7AYY7h7966yVSMyMhJPnjyBRCJBp06d4OnpCU9PTzRWc4dyUhLg4AD8f7e7CrEYSEigPmyhiY6OhpubG2QyGczMzDR2HErYakAJmxTnzp07CA8PR3h4OE6dOoX8/Hw4OTmhX79+6NmzJ9q1a0fN2zzBGMPt27cRERGB/fv34/jx48jJyUGzZs2UydvFxQUSiaTSxwoNLb4bxc+v0rsnWkYJW0AoYZMCjDHExsZi7969CA8PR3x8PKRSKbp37w5PT0989tlnqFOnDtdhkjJ4/fo1jh8/jvDwcOzbtw/Pnj2DlZUVPvvsM3h6esLDwwNVqlSp8P5v31Yk6Pv3gQYNQKPDBYwStoBQwib//fcfNm3ahM2bN+Pu3buoWbOm8oPd3d0dpqamXIeoFUlJilHQBUnI318xME/o8vPzcf78eezbtw/h4eG4efMmqlWrhkGDBsHHxwddunRRy5U3ESZK2AJCCVs/paamYuvWrdi4cSPi4uJgbm6OwYMHw8fHB59++qnefYCHhACjRulHM29SUhI2b96MjRs34vbt27C1tYW3tzeGDRuG1q1bU4liPUMJW0AoYeuPrKws7N69Gxs3bsTx48dhYGCAvn37YtiwYejTpw+MjY25DpET+jqQijGGCxcuYOPGjdi6dSuePXuGjz/+GD4+Phg6dCjq16/PdYhEC7SVsAW1vCYhXLl+/TrGjx+P2rVrY/jw4cjOzsbKlSuRmpqKXbt2YeDAgXqbrIGyFQPRRSKRCM7Ozli2bBkeP36MAwcOoHXr1pg/fz4aNmyIPn36YP/+/cgvqhQeIeVECZuQYrx79w5bt25F586d0apVK+zevRtTpkzB/fv3ERUVhdGjR8PS0pLrMHmBioEAhoaG6NOnDzZt2oS0tDT8+++/ePr0Kfr164ePPvoICxcuxLNnz7gOkwgYJWxCPpCcnIxZs2bBzs4O3t7ekEgkCAsLw8OHDzF37lxq5ixCgwYlX2E3aKDNaLhXrVo1+Pv7K8vLdu3aFT/++CPq1q0LX19fxMTEgHojSXlRwiYEir7IqKgofP7552jQoAH++OMPDB48GDdu3EBkZCQGDx5M86X/X1ISEBQEeHsr7pOSFKPBS7rCrmhNdV3Qvn17hISE4PHjx1iwYAHOnj2LDh06oG3btlizZg1ycnK4DpEIBCVsotfkcjnCw8PRoUMHdOnSBbdv38by5cvx+PFj/Pnnn/j444+5DpFXQkIUg8uWLAHCwhT3Dg7AmTOKfmqxGJBIVO/LU1Ndl9WoUQPTpk1DUlISDh06hHr16mH06NFo1KgRfv31V2RlZXEdIuE5SthEL+Xl5WHjxo1o1aoV+vfvDwMDAxw4cADXrl3DV199pfF1bYUoKUkxbUsuVywn+v59QADQqZNiNPj06cCQIYr7hATdm9JVWWKxGB4eHsrCOj179sTMmTNRv359/Pjjj0hPT+c6RMJTlLCJXnn79i3+/vtvNGnSBL6+vqhfvz5OnTqFU6dOoU+fPjR/tgRlGQneuLFiacgtWxT3+nxlXVTXwYccHBwQEhKCO3fuwNfXF4sXL4adnR2+/fZbPHr0SPtBE16jhE30wuvXr7Fo0SI0aNAAX3/9NT755BNcuXIFBw4cQKdOnbgOTxBoJHjZFdd1EBpa9PZ2dnb4/fff8eDBA0yZMgWhoaFo1KgRAgICcO/ePa3GTviLEjbRabm5uVixYgUaN26MWbNmYcCAAUhISMCWLVvQunVrrsMTFBoJXjaldR3cvl38a2vWrIn58+fjwYMH+Pnnn3Hw4EHY29vjm2++wdOnT7X3JggvUcImOkkul2Pr1q1o1qwZJkyYgB49eiAhIQGrVq1S+3KJ+oJGgpeNOorImJmZYdq0abh9+zbmzp2L9evXo1GjRpgzZw4yMjLUGzARDErYRKcwxnDkyBG0a9cO3t7ecHBwwJUrV7B+/Xo0bNiQ6/AErUkTGgleFursOjAxMUFQUBDu3r2L8ePHY/Hixfjoo4/w+++/03QwPUQJm+iM2NhYdO/eXbnsYXR0NPbv349WrVpxHZrO8POjkeCl0UTXQfXq1bF48WIkJSVhwIABmDp1Kpo2bYp169ZR2VMtKcsgQo1jpNJkMhkDwKKiorgORS89efKE+fj4MADs448/ZuHh4Uwul3MdFtFTiYmMicWMKa6nVW9iMWNJSZU/xq1bt9jAgQMZAObo6MjOnDlT+Z2SYq1dq/i/k0hU70NCFM9HRUUxAEwmk2k0DrrCJoKVm5uLX3/9Ffb29jhy5Aj+/fdfXL16Ff369eNsehYvvoUTTmmj68DBwQE7d+5ETEwMDAwM0LFjR/j5+SEtLa3yOycqKjOIUN0oYRNBOnnyJNq0aYPp06dj+PDhSExMREBAAKdrUJd3Kg/RXdrqOvjkk09w7tw5rFq1Cvv27YO9vT2WLVuGvLw89R5Ij/FpJTpK2ERQHj9+DG9vb3Tt2hVmZma4ePEi/vrrL85XzeLTt3DCD9oqIiORSDBmzBgkJibCy8sLkydPhpOTE06dOqWZA+oZPtUfoIRNBCE3NxdLly6Fg4MDIiIiEBISgtOnT6NNmzZchwaAX9/CiX6qUaMGVq5cifPnz8PY2BidO3fG8OHDkZqaynVogsan+gOUsAnvXb16Fc7Ozpg5cyZGjhyJxMRE+Pn5QSzmz68vn76FE/3Wrl07xMTEYPXq1Th48CCaN2+OjRs30nKeFcSn+gP8+cQj5AO5ubmYN28e2rVrh9evbeHr+wjPni3DokUWvBvMxadv4YSIxWKMGjUKCQkJ8PDwgK+vLwYMGICUlBSuQxMcPtUfEDH62lVpGRkZMDc3R1RUFDp37sx1ODrh2rVr8PPzw7Vr19C7dxgOHvwcIpEIjCkSIGOKPxa+zP9NSlIMMJPLCz8nFisGHFFhEcKV3bt346uvvkJubi6WLVsGHx8fWuimnG7fVnzm3L+v+AIeEPC/v+no6Gi4ublBJpPBzMxMYzHQFTbRqPJOc8rNzcX8+fPRrl075ObmIizsMg4eHAi5XMTrwVx8+hZeVlxMQaNpb9z4/PPPER8fT1fblcCLleg0OstbT1DhlKKVVmzgQ1evXmVt2rRhEomEff/99yw7O5sFBipeV1QRComEscBArb6lUiUlKWLy8lLcq6NIhiaU9/9GqMckhe3atYvVqlWLWVpasg0bNlCRITXQVuEUSthqQAm7sPJUe8rPz2dLly5lhoaGrEWLFuzChQvK57y8St6PlxcHb07gtFGJiw/HJMV7/vw58/b2ZgDYF198wV6+fMl1SIJGlc6IoJV1mtPz58/h6emJadOmYdKkSbh48SLatWun3JYGc6kfF1PQaNobv9SoUQObN2/G9u3bERERgTZt2iA2NpbrsEgpKGETjSjLNKfo6Gg4OjoiNjYWBw8exJIlSyCVSlW25dOUCl3BxRQ0mvbGT4MGDcLly5dhbW2NTp06YenSpZAXNXKS8AIlbKIRJV8ZMzx6dBpdu3ZF48aNceXKFfTu3bvIbYU4mIvvuGi1oJYS/mrQoAFOnTqFKVOmYPr06fjss8/w7NkzrsMiRaBpXWpA07oKK36aEwMgB+CAOXN8MGvWrDLV/y5pSgUpHy6moNG0N2E4dOgQhg8fDiMjI2zevBlubm5chyQINK2LCFpRV8ZisRyAHObmU3HixD/48ccfy7xYBy+mVOgILlotqKVEGHr37o2rV6+iSZMm6NatG+bOnUtN5DxiwHUARHf5+QGdOgH//stw+PB/uHp1Lzp1SsDOnYtQq1YtjRwzKUkxwKngStzfX5EsiKqC/xtttlpwcUxSfra2toiIiMD8+fMxd+5cXLx4EZs2bdLolSMpG2oSVwNqEi/e27dvMWbMGGzcuBHfffcd5s+fX6gGuLqSbEiIYsWsgkpofKyIRoiQHDx4EN7e3qhTpw7Cw8PRmL5dFYmaxIngPX78GG5ubtixYwc2b96MBQsWFErW6lpDmpa3JET9+vTpg9jYWOTl5cHZ2RnHjx/nOiS9RgmbaERsbCzat2+PJ0+e4NSpU/D29i60jTqTLM3zJUQzHBwclH/PHh4eWLZsGa38xRHBJezly5ejQYMGMDY2houLC86fP1/stl26dIFIJCp069u3r3IbPz+/Qs97eHho463orA0bNsDNzQ3169cvVAjlfepMsjTPlxDNsbS0xIEDBzBp0iRMmjQJo0ePRk5ODtdh6R1BJext27ZhypQpmDNnDi5duoTWrVujV69eePr0aZHb79q1CykpKcrbjRs3IJFIMHjwYJXtPDw8VLbbsmWLNt6OzsnPz8eMGTMwfPhweHt74+TJk7CxsSl2e3UmWZrnS4hmGRgY4JdffkFoaCg2bNiAbt26IS0tjeuw9IqgEvavv/6K0aNHY+TIkWjevDlWrlyJqlWrYu3atUVuX716ddjY2Chvx44dQ9WqVQslbKlUqrKdpaWlNt6OTsnJyYGXlxd++eUX/Pbbb1i7dm2hqmUfUmeSpYpohGjHiBEjEBUVhbt376Jjx464c+cO1yHpDcEk7Hfv3iEuLg7u7u7Kx8RiMdzd3RETE1OmfaxZswZeXl4wMTFRefzkyZOoVasW7O3tMW7cOKSnp5e4n5ycHGRkZKjc9FlGRgZ69+6Nffv2YefOnZg8eXKZ1tpVZ5Kleb6EaM8nn3yCmJgYSCQSdOjQAZcvX+Y6JL0gmIT9/Plz5Ofnw9raWuVxa2trpKamlvr68+fP48aNGxg1apTK4x4eHli/fj0iIiKwaNEiREVFoXfv3sjPzy92X8HBwTA3N1fe6tWrV7E3pQNSU1Ph5uaGS5cu4ejRoxgwYECZX6vuJOvnp6iYNX06MGSI4j4hQRhTumidaCI0DRo0wOnTp1G/fn24ubnhxIkTXIek+zS6FpgaPX78mAFgZ8+eVXl8+vTpzNnZudTXjxkzhrVs2bLU7e7cucMAsOPHjxe7TXZ2NpPJZMpbcnKyXi6vmZSUxBo2bMhsbW3ZtWvXKrEfYawhrSm0TjQRsszMTNazZ09maGjItm3bxnU4nNDW8pqCqXRmZWUFiURSaJBDWlpaiQObAOD169fYunUr5s2bV+pxGjVqBCsrK9y+fRvdu3cvchupVFpq/6yui4uLQ+/evVG9enVERkaifv36Fd5XQdlRffT+1LYPBQQoKoNRcz7hM1NTU+zbtw/+/v7w8vLC06dPMXHiRK7D0kmCaRI3MjKCk5MTIiIilI/J5XJERETA1dW1xNdu374dOTk5GDZsWKnHefToEdLT01G7du1Kx6yrjh8/ji5duqBhw4bKJjFSMTR/XFio66JoRkZGWL9+PSZPnoyvv/4aP/zwA83V1gSNXr+r2datW5lUKmWhoaEsPj6ejRkzhllYWLDU1FTGGGO+vr4sMDCw0Os6derEvvzyy0KPZ2ZmsmnTprGYmBh27949dvz4cda2bVvWpEkTlp2dXea4ZDKZ3jSJ79q1ixkaGjIPDw+WlZXFdTiC5+WlaP5WDLNTvYnFiucJP1DXRenkcjlbvHgxA8DGjh3L8vPzuQ5JK6hJvAhffvklnj17htmzZyM1NRWOjo44fPiwciDaw4cPC5W+TEhIwOnTp3H06NFC+5NIJLh27RrWrVuHV69ewdbWFj179sT8+fP1vsm7KNu3b4e3tze++OILbNy4EYaGhlyHJHg0f1wYqOuibEQiEaZPn44aNWpg1KhRyMvLwz///FPoc5lUkEa/DugJfbjC3rJlC5NIJGzo0KEsNzeX63B0RmJiyVfY+jYAj68CAxVX1EX9P0kkiueJqnXr1jGxWMz8/PxYXl4e1+FoFF1hE97YtGkThg8fjmHDhmHt2rVlXsOalK5galtAQNGrjNFVGz9Q6dvyGz58OAwMDODr64u8vDyEhobSZ0clUcImJdqyZQuGDx+OESNGYPXq1fQHpwG0TjT/UddFxQwdOhQGBgYYOnQoxGIxfeGvJErYpFhhYWEYNmwYhg0bhn///Zf6oTRIn6e2CYG/P7B4cdHPUenbkg0ZMgRyuRw+Pj4wMDDA6tWr6bOkgihhkyLt3LkTQ4cOhZeXF9auXUt/YESvUddF5Xh5eSEvLw/Dhw+HRCLBypUr6TOlAihhk0KOHj0Kb29vDBo0COvWraMmLEJAXReVNWzYMOTl5cHf3x/VqlXDL7/8wnVIgkMJm6i4cOECBg4cCHd3d2zYsAEGBvQrQkgB6rqoHD8/P2RmZuKbb75B7dq1MW3aNK5DEhT6NCZKiYmJ6NOnD1q0aIHt27fTPGtCiNp9/fXXSE1NxfTp01GrVi0MHz6c65AEgxI2AQCkpKSgV69esLKywoEDBwotQUoIIery008/ITU1Ff7+/rCyskKfPn24DkkQqNefQCaToXfv3sjNzcWRI0dQo0YNrkMihOgwkUiEVatWoU+fPhg8eDDOnTvHdUiCQAlbz2VnZ2PAgAF48OABDh8+DDs7O65DIoToAQMDA2zduhVt2rRB3759cevWLa5D4j1K2HosPz8fw4YNw7lz57Bv3z60aNGC65AIIXqkatWqCA8PR+3atdGrVy88fvyY65B4jRK2HpsxYwZ2796Nbdu2oVOnTlyHQwjRQ9WrV8fhw4cBAH369EFWVhbHEfEXJWw9FRoail9//RW///47PD09uQ6HEKLH6tatiwMHDuDOnTvw8/ODvKhl0QglbH0UExODsWPHYtSoUZg4cSLX4RABSkoCgoIAb2/FfVIS1xERoWvZsiU2btyInTt3Yv78+VyHw0s0rUvPPHr0CAMHDoSzszOWL18OUXErGhBSjJAQxdrQ75foXLxYUQHMz4/r6IiQDRgwAPPnz8esWbPQsmVLDBw4kOuQeIWusPXI27dvMWDAABgaGmLnzp0wMjLiOiQiMElJimQtlwP5+ar3AQHA7dtcR0iE7vvvv8fgwYPh6+uLq1evch0Or1DC1hOMMQQEBCA+Ph579+5FrVq1uA6JCNDatSUvM7lmjXbjURdq4ucPkUiEkJAQNG3aFP3798ezZ8+4Dok3KGHricWLF2PLli1Yt24d2rRpw3U4RKDu31c0gxeFMcXzQhMSAjg4AEuWAGFhinsHByA0lOvI9JeJiQn27t2LN2/eYNCgQXj37h3XIfECJWw9cOjQIQQFBeGHH37A4MGDuQ6HCFiDBiVfYTdooM1oKo+a+PnLzs4Ou3btQkxMDL799luuw+EFStg67tGjR/D19UWfPn0wd+5crsMhAufvX/IVdkCAduOpLE018VMTu3p06tQJy5Ytw99//43t27dzHQ7nKGHrsLy8PAwdOhTGxsYIDQ2lBeNJpTVpokhiYjEgkajer1kjvLWhNdHET03s6jV27FgMHjwYo0aNwt27d7kOh1P0Ca7D5s2bhzNnzmDLli2wsrLiOhyiI/z8gIQEYPp0YMgQxX1CgjCndKm7iZ+a2NVPJBJh9erVqFGjBry8vPS6P5sSto46ceIEfvrpJ8ydOxeffvop1+EQHdO4MRAcDGzZorgX2pV1AXU38evqKHqumZubY9u2bbhy5QqCgoK4DoczlLB1UFpaGnx8fNCtWze9/uUmpDTqbuLXxVH0fNG+fXssWrQIv/76Kw4cOMB1OJygSmcClZSk+DZ//76i2c7fX/HhI5fLMXz4cMjlcmzcuBESiYTrUAnhNT8/oFMnRYIu+HsKCKhYq4GujaLnm8mTJ+PEiRMYMWIErly5grp163IdklaJGCvu+yApq4yMDJibmyMqKgqdO3fW+PGKKg3JmOIDJzV1Ib777jscOXIEPXr00HgshJD/SUpSDDArau0KsVjR1y/U7gO+SE9Ph6OjIxo2bIgTJ07AwID7687o6Gi4ublBJpPBzMxMY8ehJnGBKXlQC8P334cgMDCQkjUhHNC1UfR8VKNGDWzZsgVnzpzBokWLuA5HqyhhC0xJg1rk8nxYWc2k+daEcKhgFH1AAGBvDzRtCoweDXTsyHVkuqNTp06YOVPxWXf9+nWuw9EaStgCU9KgFgBo2/ZzGBoaai0eQkhhp04B//4LJCYqkve//9JcbHWbM2cOmjRpAj8/P+Tm5nIdjlZQwhaYkge1iOHoaKnVeAghqmgutnZIpVKEhobi6tWrWLx4MdfhaAUlbIEpft4og0gkElxpSEJ0Dc3F1p727dtjxowZetM0TglbYD4c1CISyQHk/f+gFpFaB7XwoR4yH2IgpDxoLrZ26VPTOCVsASoY1OLrmwIgDJ06nUNCgkitpSH5UA+ZDzEQUl40F1u79KlpnOZhq4G252EDQE5ODtq1awcDAwOcP39erQPN+DCXlA8xEFIR9LvLje+++w5Lly5FXFwcWrZsqdVj0zxsUqKlS5fi1q08tG17GMOHG6q1uZgPfXB8iEEfURdE5dFcbG4UNI2PGTMG8qK+LekA7kvEkHJ78OAB5s59CLn8JtatEyurnS1erPhAqGzTOB/64PgQg74pqoKeun6n9I06y52SspFKpVixYgXc3Nywbt06jBw5kuuQ1I4StgCNGbMEubl/AxAjP1/1uYAAxQdFZT4Y+NAHx4cY9Mn7U5E+pI7fKX1UsKIZ0Z7OnTtj6NChmDlzJgYMGABLS92a5iq4JvHly5ejQYMGMDY2houLC86fP1/stqGhoRCJRCo3Y2NjlW0YY5g9ezZq166NKlWqwN3dHUk8bgc8cuQIjh6tA7G46GymjuZidS85KNQY9Al1QRBdsWTJErx9+xazZs3iOhS1E1TC3rZtG6ZMmYI5c+bg0qVLaN26NXr16oWnT58W+xozMzOkpKQobw8ePFB5fvHixVi2bBlWrlyJ2NhYmJiYoFevXsjOztb02ym3nJwcfP3116hZ0xlA0Z+u6mgu5kMfHBcx6HP/LXVBEF1ha2uLH3/8EStWrMCVK1e4Dke9mIA4OzuzCRMmKH/Oz89ntra2LDg4uMjtQ0JCmLm5ebH7k8vlzMbGhi1ZskT52KtXr5hUKmVbtmwpc1wymYwBYFFRUWV+TUX8/PPPTCKRsICAp0wiYUzxUap6k0gYCwxUz/GSkhT78vJS3CclqWe/fIxh7VrGxGLF+Xv/PiREM8fjm8BAppXfKUK04d27d6x58+asQ4cOLD8/X+PHi4qKYgCYTCbT6HEEk7BzcnKYRCJhu3fvVnl8+PDhzNPTs8jXhISEMIlEwuzs7FjdunWZp6cnu3HjhvL5O3fuMADs8uXLKq/r3Lkz++abb4qNJTs7m8lkMuUtOTlZ4wn7wYMHrGrVquzbb79liYmKZFLUh6tYzE1iFTI6n3QOiO45ceIEA8BCQ0M1fixtJWzBNIk/f/4c+fn5sLa2Vnnc2toaqampRb7G3t4ea9euxd69e7Fx40bI5XJ06NABjx49AgDl68qzTwAIDg6Gubm58lavXr3KvLUymTp1KszMzPDjjz/yoslal1D/LT+6QQhRp65du8LLywszZszAq1evuA5HLXR6lLirqytcXV2VP3fo0AHNmjXDqlWrMH/+/ArvNygoCFOmTFH+nJGRodGkHR0djR07dmDDhg3KSfk0bUR9qP9WgX6niK5ZunQp7O3tsWDBAixZsoTrcCpNMAnbysoKEokEaWlpKo+npaXBxsamTPswNDREmzZtcPv/l8speF1aWhpq166tsk9HR8di9yOVSiGVSsv5DiqGMYbAwEC0bdsWQ4cOVXmOpo2oh7qmkCUlKa7WC5Kdv7/iylVI6HeK6JI6depg6tSpWLRoESZNmoS6detyHVKlCKZJ3MjICE5OToiIiFA+JpfLERERoXIVXZL8/Hxcv35dmZwbNmwIGxsblX1mZGQgNja2zPvUtPDwcMTExGDhwoUQiwXz3yUo6phCRnXPCeGnqVOnolq1avjxxx+5DqXyNNpDrmZbt25lUqmUhYaGsvj4eDZmzBhmYWHBUlNTGWOM+fr6ssD3hrPOnTuXHTlyhN25c4fFxcUxLy8vZmxszG7evKncZuHChczCwoLt3buXXbt2jfXv3581bNiQvX37tsxxaWqUeF5eHmvevDnr3r27WvdLCgsJqfgocRqwRQi//f7770wsFrNbt25pZP/aGnQmmCZxAPjyyy/x7NkzzJ49G6mpqXB0dMThw4eVg8YePnyochX68uVLjB49GqmpqbC0tISTkxPOnj2L5s2bK7eZMWMGXr9+jTFjxuDVq1fo1KkTDh8+XKjAChc2bNiA+Ph4hNJlmsZVpv+2LIPWqJmZEO589dVX+O233/D9999j586dXIdTYbRalxpoYrWu7OxsNG3aFC4uLti+fbta9kk0w9tb0Qxe3OpMQ4YAW7ZoPy5CyP+sX78eI0aMwLlz5+Di4qLWfdNqXXpuxYoVePLkCX766SeuQyGloLrnhPCfj48PWrRogcDAQAj1OpUSNg/JZDIsWLAA/v7+sLe35zocUgqqe04I/0kkEvz88884efIkjh49ynU4FUIJm4d+++03vH79GnPmzOE6FFIGVHSEEGH47LPP0LFjR3z33XeCvMoW1KAzfZCZmYlly5Zh7NixqFOnDtfhkDKioiOE8J9IJMLcuXPh7u6Oo0ePolevXlyHVC6UsHnmn3/+QVZWFqZOncp1KKScqOgIIfzXrVs3tG/fHsHBwYJL2NQkziM5OTn49ddfMWzYMK3UJyeEEH0jEokQFBSEqKgoxMTEcB1OuVDC5pH169cjJSUFM2fO5DoUQgjRWf3790ezZs0QLLAmMUrYPJGfn4/Fixdj4MCBNDKcEEI0SCwWY+bMmdi3bx9u3LjBdThlRgmbJ3bs2IHbt28jKCiI61AIIUTnDR06FHZ2dli4cCHXoZQZJWwtSUoCgoIUVbGCghQ/F2CMITg4GD169ICTkxN3QRJCCMdK+qxUJ0NDQ0ybNg1bt27FvXv3NHMQNaOErQWlreR0+PBhXL16la6uCSF6Tdur3gUEBKB69eqCWSubEraGJSUBo0Yp6kzn56veBwQAt28Dv//+O9q3b48uXbpwHS4hhHCiLJ+V6la1alV88803CA0NxcuXL9V/ADWjhK1hpa3ktGTJCxw9ehRff/01RMVtSAgpRFtNp0Q7yrLqnSaMHj0aeXl5glgVkRK2ht2/X3Kd6aioB6hRowYGDx6s1bgIETJtN50SzSvts/L+fc0c19raGoMGDcKKFSsgL2rJPR6hhK1hJa/kxPDgwUkEBATwYv1tQoSAi6ZTonlcrno3fvx4JCUlISIiQnMHUQNK2BpW0kpOcjlDdvbfGDt2rHaDIkTAuGo6JZrF5ap3HTt2RMuWLfH3339r7iBqQAlbw4pfyYnBzm4e+vRpikaNGnEdJiGCwVXTKdEsLle9E4lEGD9+PMLDw5GcnKy5A1USJWwt8PMDEhKA6dOBIUMU92Fh1/DgwVyMHz+e6/AIERQum06JZhX1WZmQoHhc03x8fGBiYoJ//vlH8werIBET4qKgPJORkQFzc3NERUWhc+fOxW6XlKRozrt/H7hxYz9evvwFDx4ch0Qi0V6whAhcUpJigFlR44PEYsUHPC1rSipi4sSJ2LFjBx4+fAgjI6Myvy46Ohpubm6QyWQwMzPTWHx0ha0lqqNaGW7c8MCTJxHYsIGSNZ/R1CH+4bLplOi2cePGIS0tDbt37+Y6lCLRFbYalHaFXfwVAYNYLKIrAp4KCVGMRhaJFH2jBfdr1miniY6U7PZtxf/F/fuKZvCAAPo7IpXXsWNHmJub4+DBg2V+DV1h65DiR7WKaFQrT9HUIf5r3BgIDga2bFHcU7Im6jBs2DAcPXoUT58+5TqUQihhawGNahUemjpEiH4aMmQIRCIRtm3bxnUohVDC1gIa1So89CWLEP1Uo0YN9OnTB5s2beI6lEIoYWvB/woCFM4Ami4IQCqGvmQRor98fHwQGxuLJJ6NMqWErQXvj2oFciEWy2lUK89xWXWJEMKtfv36oVq1ary7yqaErSV+fsCIEcGoWnUFBg/WbkEAUn40dYgQ/VWlShUMGjQImzZtAp8mUhlwHYC+yM/Px6FDf2L06C/x++/0PUkI/PyATp1o6hAh+mjYsGEICQnB+fPn4eLiwnU4AChha01kZCRSU1Ph4+PDdSikHAqmDhFC9Iubmxvq1KmDzZs38yZh06WeluzZswf169dHu3btuA6FECIwVHFP+yQSCQYOHIg9e/bwplmcErYWMMYQHh4OT09PiIobekwIIUVQLWusuHdwAEJDuY5M9/Xr1w8PHz7EtWvXuA4FACVsrbh27RqSk5Ph6enJdSiEEAGhinvccnNzQ7Vq1bBv3z6uQwFACVsrwsPDYWZmVuJKXoQQ8iGquMctIyMj9O7dG+Hh4VyHAoAStlaEh4fDw8OjXMu1EUIIVdzjnqenJy5cuIAnT55wHQolbE17/PgxLl68SM3hhJByo4p73OvduzckEgn279/PdSiUsDWlYFTngAFvIRItRNOmfbkOiRAiMFRxj3vVq1dHp06deNEsTglbA94f1XnxYkMwNhWffGJBozoJIeVCFff4wdPTE8ePH8fr1685jUNwCXv58uVo0KABjI2N4eLigvPnzxe77erVq/Hpp5/C0tISlpaWcHd3L7S9n58fRCKRys3Dw6PC8X04qhOQADDQ6qhOmrNJiO7w81OUMZ4+HRgyhMoac8HT0xM5OTk4fvw4p3EIqtLZtm3bMGXKFKxcuRIuLi74/fff0atXLyQkJKBWrVqFtj958iS8vb3RoUMHGBsbY9GiRejZsydu3ryJOnXqKLfz8PBASEiI8mepVFrhGMsyqlOTlbNCQhRfGEQiRZOZSAQsXqw4Lv2BEyJMVHFPs5KSFJ/dBSWI/f0VrRsFGjdujIYNG+LEiRPo378/V2ECTECcnZ3ZhAkTlD/n5+czW1tbFhwcXKbX5+XlsWrVqrF169YpHxsxYgTr379/peKSyWQMAIuKimJeXoyJxYwp0qXqTSxmzMurUocqUWJiycdOStLcsQkhRIjWrlV8PkokqvchIarbjRw5krVq1arIfURFRTEATCaTaTRWwTSJv3v3DnFxcXB3d1c+JhaL4e7ujpiYmDLt482bN8jNzUX16tVVHj958iRq1aoFe3t7jBs3Dunp6SXuJycnBxkZGSq3AlyO6qQ5m4QQUnblKUzTtWtXXLt2Dc+fP+csXsEk7OfPnyM/Px/W1tYqj1tbWyM1NbVM+5g5cyZsbW1Vkr6HhwfWr1+PiIgILFq0CFFRUejduzfyFR3QRQoODoa5ubnyVq9ePeVz/xvVWXhop6ZHddKcTUIIKbvyXOS4ubkBAKKjo7UQWdEEk7Ara+HChdi6dSt2794NY2Nj5eNeXl7w9PREy5YtMWDAAOzfvx8XLlzAyZMni91XUFAQZDKZ8pacnKx8rmBUp+KXIBdiMdPaqE6as0kIIWVXnoscOzs7NGrUCJGRkdoIrUiCSdhWVlaQSCRIS0tTeTwtLQ02NjYlvnbp0qVYuHAhjh49ilatWpW4baNGjWBlZYXbJQznlkqlMDMzU7m9z88PGDv2V1StukKrozppziYhhJRdeS9yunbtWuLFnKYJJmEbGRnByckJERERysfkcjkiIiLg6upa7OsWL16M+fPn4/Dhw2Va2vLRo0dIT09H7dq1KxXv1as78dlnZ7BliwjBwUVfWat7+hXN2SSEkLIr70VOly5dcOPGDTx79kzzwRUdlHBs3bqVSaVSFhoayuLj49mYMWOYhYUFS01NZYwx5uvrywIDA5XbL1y4kBkZGbEdO3awlJQU5S0zM5MxxlhmZiabNm0ai4mJYffu3WPHjx9nbdu2ZU2aNGHZ2dlljuv9UeIF+zUwMGB///13sa8p68jEikhKYiwwUDEiPTCQRocTQkhxQkLK/lmcnJzMALDt27erPK6tUeKCStiMMfbnn38yOzs7ZmRkxJydndm5c+eUz7m5ubERI0Yof65fvz6DYvSXym3OnDmMMcbevHnDevbsyWrWrMkMDQ1Z/fr12ejRo5VfAMrqw4R95MgRBoDFx8cXuT1NvyKEEP4oz0VO48aNVaYXM6a9hC2owikAMHHiREycOLHI5z7sW7hfyrDoKlWq4MiRI2qK7H9iY2NhYWEBBweHIp/nurgKIYSQ/ylPYZoOHTogNjZWswEVQzB92EISFxcHJycniIrJyjT9SrdQKVhC9IeTkxOuXbuGd+/eaf3YlLA1IC4ursQBbjT9Sne8v9BLWJji3sEBtNALITrKyckJ7969w82bN7V+bErYavb06VM8evQITk5OxW5D0690Q3mqJJGyodYKwneOjo4Qi8WIi4vT+rEpYatZwX9iSQmbpl/pBioFq17UWkGEwMTEBA4ODpwkbMENOuO7uLg4WFpaomHDhiVu5+cHdOqk+FAvWCEmIICStZDQWAT1eb+14kMBAYq/FfrbIHzRrl07Sti64OLFi2jbtm2xA87eR0vmCRuNRVAfmjlBhMTJyQlbt27Fu3fvYGRkpLXjUpO4mpU24IzoDhqLoD7UWkGEhKuBZ5Sw1ejly5elDjgjuoPGIqgPtVYQIeFq4Bk1iavR3bt3AQCtW7fmOBKiLTQWQT38/YHFi4t+jlorCN+YmJigcePGuH79ulaPSwlbjZKTk2FgYFDqgDOiW2gsQuUVtFYEBCiuqBn73z21VhA+sre3R2JiolaPSQlbjR49eoRGjRrB0NCQ61AIEZzytlYkJSkGqxVs6++vSPyEaIO9vT127dql1WNSwlajhw8fwt7enuswCBGssrZWhIQopoG9fzW+eLEi2Wt63XlCAKBp06a4f/8+cnJytHZMGnSmRo8ePULTpk25DoMQnUYV5ggfNG3aFHK5XDl2SRsoYatRSkoKJWxCNIwqzBE+KGhNTUhI0NoxKWGrkVwup4RNeEtX6nTTnG3CB9bW1qhWrZpWB55RwlYzStiEj3SpTjfN2SZ8IBKJ0LRpU0rYQlWlShXUrl2b6zAIUaFrfb5UYY7whbandlHCVqM6deqUqYY4Idqka32+VGGO8EWTJk2QpMW+JZrWpUY1a9bkOgRCCtHFPl+qMEf4oE6dOkhLS0NeXp5WjkcJW42srKy4DoGQQnS1z5cqzBGu1a5dG4wxvHr1SivHoyZxNapRowbXIRBSCPX5EqIZtra2AID09HStHI8SthpZWVnpzNQZojuoz5cQzSgYZKythE1N4mp0714HODhQuUTCP9TnS4j61apVC2KxWGsJW8RYcY1lpKwyMjJgbm4O4BUA80LPi8VAQgJ9OBJCiK6xtbWFh4cHQkJCIJPJYGZmprFj0RW2GolErMi+woKpM+oYICPkFYqEHDshhBSldu3aeP78uVaORQlbrYoeiquuqTNCXqFIyLETQkhxbG1tadCZEGly6oyQq1UJOXZCCClJ7dq1KWELkSanzgi5WpWQYyeEkJLUqFEDGRkZWjkWJWw1mjnzocamzgi5WpWQYyeEkJKYmpri7du3WjkW9WGrUe/eLxEQoJmpM0KuViXk2AkhpCTVqlXDmzdvtHIsSthqpqlyif7+ikFaReF7tSohx04IISWpVq0a3r17p5VjUZO4QAi5WpWQYyeEkJKYmppq7Vh0hS0gQq5WJeTYCSGkONWqVdPasShhC4w2VyhSd6ETWl2JEKJrKGGTctFEBTEqdEIIIaWjJnFSZppIrO8XOvlQQICiaZuasglRHyrbK1zavMKmQWcCpqkKYlTohBDtCQkBHByAJUuAsDDFvYMDEBrKdWRFoyWEVZmYmGjtWJSwBUxTiZUKnRCieUlJwNixiqtpoZTtFdqXC20wMNBeQ7XgEvby5cvRoEEDGBsbw8XFBefPny9x++3bt8PBwQHGxsZo2bIlDh48qPI8YwyzZ89G7dq1UaVKFbi7uyNJIF8ZNZVYqdAJIZpVkPhWry5+G761ZtGaAEUTi7WXRst9pBEjRiA6OloTsZRq27ZtmDJlCubMmYNLly6hdevW6NWrF54+fVrk9mfPnoW3tzcCAgJw+fJlDBgwAAMGDMCNGzeU2yxevBjLli3DypUrERsbCxMTE/Tq1QvZ2dnaelsVpqnE6u+v2brohOiz9xNfcX9nAP9as6irrGjaTNhg5dS/f39maGjIGjduzBYsWMAePXpU3l1UmLOzM5swYYLy5/z8fGZra8uCg4OL3H7IkCGsb9++Ko+5uLiwsWPHMsYYk8vlzMbGhi1ZskT5/KtXr5hUKmVbtmwpNo7s7Gwmk8mUt+TkZAaARUVFVebtlVtiImNiMWOKP23Vm1jMWFJSxfcdEqLYh0Sieh8Soq7oCdFPgYGKv6ei/m7fv0kkim35wsur5M8bLy+uI+RGRkYGA8AAMJlMptFjlfurwZ49e/D48WOMGzcO27ZtQ4MGDdC7d2/s2LEDubm5av9CUeDdu3eIi4uDu7u78jGxWAx3d3fExMQU+ZqYmBiV7QGgV69eyu3v3buH1NRUlW3Mzc3h4uJS7D4BIDg4GObm5spbvXr1KvPWKkyTFcT8/ICEBGD6dGDIEMV9QgJN6SKkskrqynof31qzqKusaKLiTooGVOhavmbNmpgyZQquXr2K2NhYNG7cGL6+vrC1tcW3336rkT7g58+fIz8/H9bW1iqPW1tbIzU1tcjXpKamlrh9wX159gkAQUFBkMlkyltycnK534+6aDKxFhQ62bJFcU9TuQipvJISH6B4jo9leyvTVabLI8vlRc1/1ZBKNb6npKTg2LFjOHbsGCQSCfr06YPr16+jefPm+O2339QVI+9IpVKYmZmp3LhEiZUQ4Sgp8QHAmDH8bM2qaIuero8s53XCzs3Nxc6dO/HZZ5+hfv362L59OyZPnownT55g3bp1OH78OMLCwjBv3jy1BmplZQWJRIK0tDSVx9PS0mBjY1Pka2xsbErcvuC+PPskhJDKKCnxhYQAK1fy90t3eVv09GFkuTYTdrknkNWuXRtyuRze3t44f/48HB0dC23TtWtXWFhYqCG8/zEyMoKTkxMiIiIwYMAAAIoTFRERgYkTJxb5GldXV0RERGDy5MnKx44dOwZXV1cAQMOGDWFjY4OIiAjl+8jIyEBsbCzGjRun1vgJ4QOqqMUPQl4MpzxrApRlZLnQ1xfIz8/X3sHKO0pt/fr17O3bt5oYAFeqrVu3MqlUykJDQ1l8fDwbM2YMs7CwYKmpqYwxxnx9fVnge8Mqz5w5wwwMDNjSpUvZrVu32Jw5c5ihoSG7fv26cpuFCxcyCwsLtnfvXnbt2jXWv39/1rBhw3K9R5lMxskocaI+iYmKEbleXor7xESuI1K/tWtp5D/RLn0YWV4wSwhaGCVe7oTNtT///JPZ2dkxIyMj5uzszM6dO6d8zs3NjY0YMUJl+7CwMNa0aVNmZGTEPv74Y3bgwAGV5+VyOZs1axaztrZmUqmUde/enSUkJJQrJkrYwqYPiUyTUwAJKU5JU9j4Nm2touLj47WWsEWMlWWCASlJRkYGzM3NERUVhc6dO3MdDimHpCTFAJiiuqHEYkX/nBCaKUsTFKQY7FNU651EouiLFHrTJOEfffj7On/+PFxcXAAAMplMo4OQBVealBB10pfqTVQfnnBBk7Ui+CIrK0trx6LlNYnW8Wngk74kMip6Qbgi5AF2ZZGZmam1Y1HCJlqlifW7K0NfEpm/v+I8F4VvFbWI7inPyHKh0WbCpiZxojV8nJOp6YVO+FLhSR+aJgnhQlZWltbKk1LCJlrDx/5iTSYyvlV4ovrwhKhfZmYmqlatqpVjUZM40Rq+9hdroo/t/daEDwUEKI7HxVWtLjdNEsKFgoT9+vVrjR+LEjbRGj73F6s7kelDhSdCiGJar4mJCZ49e6bxY1GTuBrRlPaSabq/mE/42ppACFGv1NRUVK9eXSvHooStRtpoEhEyfRr4xOfWBEKI+qSkpFDCFqL09HSuQ+A9fRn4pE+tCYTosydPnqBGjRpaORb1YavR8+fPuQ5BEPRh4FNBa0JAgOqcc8Z0rzWBK3wqwEP0E2MMKSkpsLKy0srxKGGrEV1hk/fpeoUnLvGtAA/RT5mZmXj9+jVdYQsRJWzyIX1oTdA2vk6ZI/onJSUFALSWsKkPW41evHjBdQikgvhSkYyUjo8FeIh+KkjY1CQuQNSHLUzUvCosNGWO8MWTJ08A0BW2IFHCFh4+1jcnJaMpc4QvUlJSYGpqiipVqmjleJSw1aigeYQIBzWvCg9NmSN8cffuXdSvX19rx6OErUbPnj2j4ikCQ82rwqNPBXgIvyUmJsLe3l5rx6OErWa3qQ1VUKh5VZj0pQAP4bfExEQ0bdpUa8ejQWdqlpiYiNatW3MdBikjf3/FALOiaLN5lYqAlB9NmSNcevv2LR4+fKjVhE1X2GpkZmaGxMRErsMg5cCH5lW+rZtNCCld0v/P/dRmkzhdYatRvXr1KGELEJcVyagICCHCVPBZ37RpU8THx2vlmJSw1ahu3bpISEjgOgxSAVw1r+riutnUvE/0QWJiIiwtLbU2BxughK1WdnZ2iI2NBWMMouI+hQl5j66NUq9sERpK9kQoEhIS0LRpU61+1lMfthrVq1cPr169opripMx0aZR6ZYvQUF8+ERJtT+kCKGGrlZ2dHQDg5s2bHEdChEKXioBUpggNVZwjQsIYw61btyhhC1m9evVQpUoVxMXFcR0KEQg+jFJXl8o071PFubKhRWr44c6dO5DJZGjbtq1Wj0t92GpkYGAAR0dHnUnY1J+oHbqybnZlmvd1rS+/okr6m6NFavjj4sWLAAAnJyetHpcStpo5OTnh2LFjXIdRorIkYvpw0C5dKAJSmSI0utSXX1El/c117EjT//gkLi4OdnZ2qFmzplaPS03iaubk5ITExERkZmZyHUqRyjKwh/oTSUVUpnlfl/ryK6K0v7lffqEuAz6Ji4vT+tU1QAlb7ZycnMAYw+XLl7kOpZCyJmLqTyQVVdEa37rUl18Rpf3NRUdTlwFfMMZw6dIlThI2NYmrWbNmzVClShVcvHgRnTt35jocFWUt0kH9iaQyKtq8ryt9+RVR2t8cQF0GfFEw4IwStg7g88CzsiZi6k8kXNGFvvyKKO1v7tNPFS0VRdGHLgM+KfhspyZxHeHk5FShhK3pKRtlTcT63p9IiLaV9jc3fbp+dxnwCVcDzgBK2BrRrl07JCYm4tWrV2V+jTaqPJU1Eet7fyIh2laWvzlaA5wfLly4wMnVNUBN4hrx6aefgjGG06dP47PPPit1e22t2FTwoRAQoDp1hLHCiVif+xMJ4UJZ/ub0tcugrDRdOyI7OxsxMTFYuHCh+nZaDpSwNaBhw4aoV68eIiMjy5SwtbliU3kSMX04EKJd9DdXcdqoHREbG4ucnBx06dJFPTssJ8E0ib948QI+Pj4wMzODhYUFAgICkJWVVeL2X3/9Nezt7VGlShXY2dnhm2++gUwmU9lOJBIVum3durVSsYpEInTp0gUnT54s0/baHpVd8KGwZYvinq6aCV9Q6U1SEdqqHXHy5ElYWlqiVatW6tlhOQnmCtvHxwcpKSk4duwYcnNzMXLkSIwZMwabN28ucvsnT57gyZMnWLp0KZo3b44HDx7gq6++wpMnT7Bjxw6VbUNCQuDh4aH82cLCotLxdunSBRs3bsSrV69K3Z++jMqmUqekJFRdj1SUtlopIyMj4ebmBrGYo2tdJgDx8fEMALtw4YLysUOHDjGRSMQeP35c5v2EhYUxIyMjlpubq3wMANu9e3e54snOzmYymUx5S05OZgBYVFSUcpu7d+8yAGzv3r2l7i8xkTGxmDHFx5TqTSxmLCmpXOHx0tq1ivcikajeh4RwHRnhA334GyCa4+VV8u+Pl1flj/H27VsmlUrZ77//Xui5qKgoBoDJZLLKH6gEgmgSj4mJgYWFBdq1a6d8zN3dHWKxGLGxsWXej0wmg5mZGQwMVBsWJkyYACsrKzg7O2Pt2rVgxbVP/7/g4GCYm5srb/Xq1Su0TYMGDWBnZ1emZnE+jcrWRJMklTolpaHqeqQytNFKGRMTg5ycHHTt2rXyO6sgQSTs1NRU1KpVS+UxAwMDVK9eHampqWXax/PnzzF//nyMGTNG5fF58+YhLCwMx44dwxdffIHx48fjzz//LHFfQUFBkMlkyltycnKhbcrbj82HKRuamlpGH8akNFRdj1SGNmpHnDx5EtWrV0eLFi0qv7MK4rQPOzAwEIsWLSpxm1u3blX6OBkZGejbty+aN2+OH3/8UeW5WbNmKf/dpk0bvH79GkuWLME333xT7P6kUimkUmmpx+3SpQs2bNiAFy9eoHr16qVuz+UIUU1OLaMPY1IafRnHQTSjPFNWK4rz/mtwfIU9depU3Lp1q8Rbo0aNYGNjg6dPn6q8Ni8vDy9evICNjU2Jx8jMzISHhweqVauG3bt3w9DQsMTtXVxc8OjRI+Tk5FT6/XXt2hWMMURFRVV6X5qmyatg+jAmpaHqeqSyNNlK+fr1a8TGxnI2nasAp1fYNWvWLFN5N1dXV7x69UplSbMTJ05ALpfDxcWl2NdlZGSgV69ekEqlCA8Ph7GxcanHunLlCiwtLct0BV2aBg0awMHBAfv378fnn39e6f1pkiavgiuzTjLRD9q4QiK6T1OtlMePH8e7d+/Qu3dv9e+8HAQxratZs2bw8PDA6NGjsXLlSuTm5mLixInw8vKCra0tAODx48fo3r071q9fD2dnZ2RkZKBnz5548+YNNm7ciIyMDGRkZABQfFGQSCTYt28f0tLS8Mknn8DY2BjHjh3Dzz//jGnTpqkt9n79+mHdunWQy+WcNqWURpNXwfRhTD5U1BQ/qq5H+Co8PBwODg5owvU8VI2OQVej9PR05u3tzUxNTZmZmRkbOXIky8zMVD5/7949BoBFRkYyxhiLjIxkAIq83bt3jzGmmBrm6OjITE1NmYmJCWvdujVbuXIly8/PL1dsMpms0LSuAqdOnWIAWExMTIXfuzZoY1pNUhJjgYGKKRaBgTRVR1/RFD8iJHl5eaxmzZpsxowZxW6jrWldIsZKmcNESpWRkQFzc3NERUUVWgM7Pz8f1tbWGDNmDH7++WeOIiyb0NDir4KpcAVRh6QkxcyDogY3isWKPke6oiZ8EhMTgw4dOuD06dPo2LFjkdtER0fDzc1NOXVYUwTRJC5kEokEn332GcLDw3mfsKlJkmhaaYMblywBqlenaniEP8LDw2FlZYVPPvmE61AoYWuDp6cn1q1bhzt37uCjjz7iOpwS0eIDRJNKGtwolwP//KMoGkSlSQlf7Nu3D5999hkkEgnXoQijcIrQ9ezZE0ZGRti3bx/XoRDCqZIGNxYkcqqGR/jizp07uHnzJjw9PbkOBQAlbK0wNTVF9+7dER4eznUohHCqpPnWxaFqeIQr+/btg1QqRY8ePbgOBQAlbK3x9PREdHQ00tPTuQ6FEM4UVzcfKPnKm6rhES7s2bMH3bp1g6mpKdehAKCErTUFhVO2b9/OcSSaQ2sZk7IoqiLV2LGK5F0UqoZHuJCcnIzo6GgMGTKE61CUaNCZllhbW6NHjx7YtGkTvvrqK67DUTtay5iUx4eDG5OSgNWri96WquERLmzZsgVSqRQDBw7kOhQlusLWomHDhuH06dO4d+8e16GoFS2fSSqLT0vMEgIAmzZtQv/+/TU6r7q8KGFr0YABA2BiYoLNmzdzHYpa0fKZ6qWvXQt8WGKWEAC4du0arl27Bh8fH65DUUFN4lpkYmKCAQMGYOPGjfjuu+8gKi7LCQwtn6k++t61QHUAyPuKqjmvjUI6mzZtQo0aNdCrVy/NH6wc6Apby4YNG4b//vsPly9f5joUtaHlM9WDuhYI+Z+QEEUZ2yVLgLAwxb2Dg6KEsibJ5XJs3rwZQ4YMgZGRkWYPVk6UsLXM3d0dtWrVwqZNm7gORW1oLWP1oK4FQhS4/PIaHR2NR48eYdiwYZo7SAVRwtYyAwMDeHl5YcuWLcjPz+c6HLWgAUPqQV0LhChw+eV148aNaNiwIVxdXTV3kAqihM2BYcOGISUlBceOHeM6FLWhAUOVR10LhChw9eX19evX2L59O3x8fHg5xogGnXGgXbt2cHR0xIoVK+Dh4cF1OGpDA4Yqx99fMcCsKNS1QPQJV19et2zZgszMTATw9I+NrrA5IBKJMH78eOzfvx8PHjzgOhzCE9S1IBz6OvVOW7gYF8MYw/Lly9G3b1804GlzFiVsjgwdOhSmpqZYtWoV16EQHqGuBf7javSyPuHiy+u5c+dw5coVTJgwQf07VxMRY+VdO4d8KCMjA+bm5oiKikLnzp3L/LpJkyZhy5YtSE5OhlQq1WCERKi4modKipaUpEjOcnnh58RixZcraglRn9u3FQm64Pc/IEBz59fX1xdnz55FUlISxMUVti9GdHQ03NzcIJPJNFoZja6wOTRu3Dg8e/YMO3fu5DoUwkN0JVd56m66pql32lUwLmbLFsW9ppL1s2fPEBYWhnHjxpU7WWsTfyPTAw4ODujWrRv+/vtvrkMhPENFVCpPE194aOqdsBX3BW7t2rUQiUQYOXIktwGWghI2x8aPH48zZ87g6tWrXIdCeISu5CpHU194aOqdcBX3BW7tWjlWrlwJLy8v1KhRg+swS0QJm2Oenp6wtbWlq2yigq7kKkdTX3ioqp8wlfQFbtQoEe7fl2D8+PFch1kqStgcMzQ0xFdffYUNGzbg2bNnXIdDeIKu5CpHU194aOqdMJX0BY4xOWxsfkD79u21G1QFUMLmgQkTJkAsFuOPP/7gOhTCE3QlVzma/MJDU++Ep6QvcADDRx9152Vlsw9RwuaB6tWrY+zYsfjrr7+QkZHBdTiEB+hKrnI0/YVHW6OXiXqU9AUOADp1qqu1WCqDEjZPTJkyBW/evMHKlSu5DoXwBF3JVRx94SHvK/4LHINIJMaoUfy/ugaoljhv1KlTByNGjMBvv/2Gb775BsbGxlyHREqhjaImVJ+94vz8gE6dtFd4g/BXwRe4gADFlTZjAGP5YAz45x+Gxo2Fce1Klc7UoKKVzj6UlJQEBwcHLF++HF999ZUaIyTqFhKiGHVa8MdfcL9mDV0BE8JXBZXTbt7Mwv79y/Hdd9b46Se/Su+XKp3poSZNmmDQoEFYvHgx8vLyuA6HFIOKmhAiTAUtVnZ2gahefQmCggZzHVK5UMLmmcDAQNy7dw9hYWGV3hetKKQZVNSEEOFKS0vDmjVrMGnSJJiYmHAdTrlQwuaZNm3awMPDAz///DPy8/MrvB+qQ605VNSEEOH69ddfYWBggIkTJ3IdSrlRwuahH3/8ETdv3sTmzZsr9HpqstUsKmpCiDA9fvwYy5Ytw6RJk2Bpacl1OOVGCZuHXFxc8Pnnn2PWrFnIyckp9+upyVazqKgJIcI0b948mJiYYPr06VyHUiGUsHlqwYIFSE5OrtC8bGqy1Sya40uI8CQkJGDNmjX47rvvYG5uznU4FUIJm6eaNWsGPz8//PTTT8jMzCzXa6nJVvOoqAkhwjJr1izY2toKYpGP4lDC5rEff/wRmZmZ+OWXX8r1uso02dLI8rKj8pSECMPFixexfft2zJ07V9BFqShh81i9evUwceJE/PLLL3j69GmZX1fRJlsaWU4I0UVBQUFo1qwZfH19uQ6lUgSTsF+8eAEfHx+YmZnBwsICAQEByMrKKvE1Xbp0gUgkUrl9WEHs4cOH6Nu3L6pWrYpatWph+vTpvCpaEhQUBLFYjAULFpTrdeVtsqWR5YQQXXT8+HEcP34cP//8MwwMhF2NWzDR+/j4ICUlBceOHUNubi5GjhyJMWPGlDr1afTo0Zg3b57y56pVqyr/nZ+fj759+8LGxgZnz55FSkoKhg8fDkNDQ/z8888aey/lUaNGDcyYMQNz587FpEmT0KhRozK/tjx1qMsyslwfalproz44IUQ75HI5AgMD8cknn6B///5ch1N5TADi4+MZAHbhwgXlY4cOHWIikYg9fvy42Ne5ubmxSZMmFfv8wYMHmVgsZqmpqcrHVqxYwczMzFhOTk6xr8vOzmYymUx5S05OZgBYVFRU+d5YGWVlZbG6deuy/v37a2T/jDHm5cWYWKwoif/hTSxWPK/r1q5VvFeJRPU+JITryAghFfHvv/8yAOzUqVMaPU5UVBQDwGQymUaPI4gm8ZiYGFhYWKBdu3bKx9zd3SEWixEbG1viazdt2gQrKyu0aNECQUFBePPmjcp+W7ZsCWtra+VjvXr1QkZGBm7evFnsPoODg2Fubq681atXrxLvrnQmJib45ZdfsHfvXhw6dEgjx9D3keVC7xKgwYKEqHrx4gUCAwMxbNgwdOrUietw1EIQCTs1NRW1atVSeczAwADVq1dHampqsa8bOnQoNm7ciMjISAQFBWHDhg0YNmyYyn7fT9YAlD+XtN+goCDIZDLlLTk5uSJvq1wGDx6Mbt264euvv0Z2drba96/vxUCEXGyGBgsSUlhB4anFixdzHYracJqwAwMDCw0K+/D233//VXj/Y8aMQa9evdCyZUv4+Phg/fr12L17N+7cuVOpuKVSKczMzFRumiYSifDXX3/hwYMHWLp0qdr3r+/FQIRabEboLQOEaMKlS5ewYsUKzJ07F7Vr1+Y6HLXhdNDZ1KlT4VdKpYlGjRrBxsam0LSmvLw8vHjxAjY2NmU+nouLCwDg9u3b+Oijj2BjY4Pz58+rbJOWlgYA5dqvtjRr1gyTJ0/Gzz//DF9fX9SvX1+t+/fzAzp1UiTogkFXAQG6n6wB4XYJ0GBBQlTJ5XJMmDABzZs3F+QCHyXhNGHXrFkTNWvWLHU7V1dXvHr1CnFxcXBycgIAnDhxAnK5XJmEy+LKlSsAoPzG5erqigULFuDp06fKJvdjx47BzMwMzZs3L+e70Y7Zs2dj8+bN+Pbbb7Fr1y617788I8t1ib8/UFzLGZ+7BITaMkCIpqxbtw7nzp3DyZMnYWhoyHU4aiWIPuxmzZrBw8MDo0ePxvnz53HmzBlMnDgRXl5esLW1BaBYhcXBwUF5xXznzh3Mnz8fcXFxuH//PsLDwzF8+HB07twZrVq1AgD07NkTzZs3h6+vL65evYojR47ghx9+wIQJEyCVSjl7vyWpVq0afvnlF+zevRtHjhzhOhydIdQuAaG2DBCiCS9fvsTMmTMxdOhQuLm5cR2O+ml0DLoapaenM29vb2ZqasrMzMzYyJEjWWZmpvL5e/fuMQAsMjKSMcbYw4cPWefOnVn16tWZVCpljRs3ZtOnTy807P7+/fusd+/erEqVKszKyopNnTqV5ebmlis2mUym0WldH5LL5axr166sSZMm7O3bt1o5pr5ISmIsMFAxjS0wUPEznyUmljwdj+/xE6JOEyZMYKampiVO99UEbU3rEjFWXIMaKauMjAyYm5sjKioKnTt31sox4+Pj0aZNG0yZMgXB+tiGTZRCQxVN9iKRIlUX3K9ZQ4uREP1x6tQpuLm54ZdffsG3336r1WNHR0fDzc0NMplMo4OQBVPpjKhq3rw55syZg1mzZuHzzz+Hs7Mz1yERjujzYEFCAODNmzfw9/eHq6srvvnmG67D0RhK2AI2Y8YM7Nq1C35+frh06ZKgV6EhlaOvgwUJAYDvvvsOjx49woEDByCRSLgOR2MEMeiMFM3AwAChoaG4c+cO5s6dW+h5qn5FCNF1p06dwrJly7BgwQI0bdqU63A0iq6wBa5FixZFNo2HhCgKarzfr7l4MfVrfogW+yBEuAqawjt06IBJkyZxHY7G0aAzNeBi0Nn78vLy4OrqitevX+PSpUtITjaGg4Oi2tWHxGLFMpvUv1n0lxouB2vRlwf+of+Twvh0TiZPnoxVq1bh6tWrnF5da2vQmWCmdfGZtqd1FeX69evMyMiIBQYGssBAxUpTRU31kUgU05X0Hd+mQ9FKYfxD/yeF8emcREdHM5FIxH799VftH/wDtFoXKZeCpvHFixcjLu45Vb8qBZ8W+6B64PxD/yeF8emcZGVlKZvCdXlU+IcoYeuQGTNmwNnZGbGxYRCJis7YVP1KgU8lPfn05YEo0P9JYXw5J4wxjBs3DikpKQgJCdHpUeEfooStQwwMDLB161Ywtgb5+QxA4YzE57rY2sSnkp58+vJAFOj/pDC+nJN169Zh48aNWLVqFZro2YACStg6pn79+ggN/R6AP0QiJqi62NrEp/W/+fTlgSjQ/0lhfDgnt27dwoQJEzBy5Ej4+Pho/oA8QwlbBw0cOBATJ1aDoWELDB+ehiFDgOnTFaPDaUqXAp8W++DTlweiQP8nhXF9Tt6+fYsvv/wS9evXx59//qnZg/EUzcPWUUuWLMHp0644c6Yz4uLiYGpqynVIvMOXkp4FXx6KqwdOLSLaR/8nhXF9Tr799lskJSXh/PnzMDEx0ezBeIrmYasB1/Owi5OYmAgnJycMGDAA69evh6i49iyO8WleJ5du3+b+ywNRRf8nhXFxTsLCwvDll19i1apVGDNmjGYPVgHamodNCVsN+JqwAWDTpk0YNmwYQkJC4MfD9nC+FS8hhPDL3bt30aZNG3h4eGDr1q28vPDQVsKmPmwd5+PjA39/f4wfPx6XL1/mOhwVfJrXqQ1U252Q8nn9+jW++OILWFlZ4Z9//uFlstYmSth64K+//kLz5s3Rv39/pKWlcR2OEl/mdWpDSAjg4AAsWQKEhSnuHRwUa1kTQgpjjGHkyJFITEzErl27YG5uznVInKOErQeqVKmCPXv24N27d/jiiy/w7t07rkMCwJ95nZqmby0JhKjDggULsH37dqxfvx6tW7fmOhxeoIStJ+rWrYvdu3fjwoULmDBhAvgwdIEP8zq1QZ9aEghRhz179mDWrFn48ccf8cUXX3AdDm9QwtYjrq6u+Oeff/Dvv/9i+fLlXIfD+bxObdGXlgRC1OH69evw9fXFF198gVmzZnEdDq9QwtYzI0aMwJQpUzB58mRERERwGgufipdokr60JBBSWc+fP0f//v3x0UcfYd26dRCLKUW9j86GHlq0aBG6d++OwYMH486dO5zG4uenqMA2fTp0tiKbvrQkEFIZubm5GDx4MDIzM7F37169LY5SEqp0pocKFglxcXFBv379cObMGVhaWnIWT+PGQHAwZ4fXOK4rRBHCd4wxTJw4EadPn0ZERATq16/PdUi8RFfYesrS0hL79u3D06dP4enpibdv33Idkk7Th5YEQipq/vz5+Oeff7Bq1SreFZ/iE7rC1mP29vY4cOAAunXrBi8vL+zcuRMGBvQroSm63pJASEWsWrUKc+bMwYIFC+Dv7891OLxGV9h6zsXFBTt27MCBAwcwbtw4Xkz34guqTEaIZu3atQvjx4/H119/jaCgIK7D4T26nCLo3bs31q5dixEjRsDGxgbz58/nOiTOFVXjfPFiqnFOiLpERUVh6NChGDRoEH7//Xe9LztaFpSwCQBg+PDhSEtLw4wZM2BtbY2JEydyHRJn3q9M9qGAAMWSnDRQjJCKu3r1Kjw9PdGpUyesX7+epm+VESVsHuJquclp06YhNTUV33zzDWrVqoUhQ4Zo/qA8VJbKZNQXTUjF3L9/Hx4eHvjoo4+wa9cuSKVSrkMSDErYPMNlU6xIJMKSJUuQlpYGX19fWFpaokePHpo9KA9RZTJCNCM1NRW9evWCiYkJDh06pNGlKHURtUPwCB8WiRCLxVi7di3c3d3h6emJEydOaP6gPEOVyQhRv6dPn6J79+7IysrCkSNHYG1tzXVIgkMJm0f4skiEkZERdu7cic6dO6Nfv36IiorSzoF5giqTEaJez58/h7u7O168eIETJ07go48+4jokQaKEzSN8aoo1NjbGnj174Orqir59++L06dPaOzjH9KXGOSHa8OLFC7i7uyMtLQ0nTpyAvb091yEJFiVsHuFbU2yVKlUQHh6O9u3bw8PDA9HR0doNgENUmYyQynv+/Dm6deuGx48fIyIiAs2aNeM6JEGjQWc84u+vGGBWFK6aYqtWrYr9+/fD09MTvXv3xv79+9G1a1ftB1IMTY6op8pkhFTc06dP4e7ujtTUVERGRqJFixZchyR4dIXNI3xtijUxMcG+ffvQsWNH9O3bF8ePH+cmkA+EhAAODsCSJUBYmOLewQEIDeU6MkL0W2pqKrp27YqnT5/i5MmTlKzVhBI2z/C1KbZq1aoIDw+Hm5sb+vXrh/DwcE7j4cOIekJIYQ8ePICbmxtevnyJqKgoNG/enOuQdAYlbB4qaIrdskVxz5dBTgUD0fr27YvPP/8ca7Q1bL0IfBlRTwj5n+vXr6NDhw7Izc1FdHQ0DTBTM8Ek7BcvXsDHxwdmZmawsLBAQEAAsrKyit3+/v37EIlERd62b9+u3K6o57du3aqNtyRIUqkU27Ztw9ixYzFq1CgsWLCAkwVD+DSinhACnDp1Cp9++ilq1aqFs2fPojFfrjR0iGAGnfn4+CAlJQXHjh1Dbm4uRo4ciTFjxmDz5s1Fbl+vXj2kpKSoPPbPP/9gyZIl6N27t8rjISEh8PDwUP5sYWGh9vh1iUQiwfLly1G7dm388MMPSE1NxR9//KHVesB8G1FPiD7bs2cPvLy80KFDB+zZs4cqmGkKE4D4+HgGgF24cEH52KFDh5hIJGKPHz8u834cHR2Zv7+/ymMA2O7duysVn0wmYwBYVFRUpfYjRKtWrWJisZgNGTKEZWdna+24iYmMicWMKa6nVW9iMWNJSVoLhRC99s8//zCxWMwGDRqk1c8APomKimIAmEwm0+hxBNEkHhMTAwsLC7Rr1075mLu7O8RiMWJjY8u0j7i4OFy5cgUBRcyNmjBhAqysrODs7Iy1a9eW2sSbk5ODjIwMlZu++HCN6K5dx2DHjh3Yu3cv+vTpo7VzwdcR9YToC8YY5s+fjzFjxmDcuHHYunUrLeShYYJoEk9NTUWtWrVUHjMwMED16tWRmppapn2sWbMGzZo1Q4cOHVQenzdvHrp164aqVavi6NGjGD9+PLKysvDNN98Uu6/g4GDMnTu3/G9E4IpfmORzHD16FJ6enujSpQv2798PW1tbjcfj56dY6nLNmv/Nww4IoGRNiKbl5eVh0qRJ+PvvvzF//nx8//33tJ61Nmj0+r0UM2fOZABKvN26dYstWLCANW3atNDra9asyf7+++9Sj/PmzRtmbm7Oli5dWuq2s2bNYnXr1i1xm+zsbCaTyZS35ORknW8SL0sT9LVr11idOnVY7dq1WWxsLNchE0I0ID09nbm7uzMDAwO2evVqrsPhBW01iXN6hT116lT4lTLBuFGjRrCxscHTp09VHs/Ly8OLFy9gY2NT6nF27NiBN2/eYPjw4aVu6+Ligvnz5yMnJ6fY5h2pVKp3TT9lWyO6JS5evIiBAweic+fO+PfffzFs2DDtBkoI0Zj4+Hh4enri1atXOHbsGLp06cJ1SHqF04Rds2ZN1KxZs9TtXF1d8erVK8TFxcHJyQkAcOLECcjlcri4uJT6+jVr1sDT07NMx7py5QosLS31LiGXpqzTqGxsbBAZGYmvvvoKvr6+uHbtGoKDgyGRSLQVKiFEA/bt2wcfHx/Ur18fx44dQ8OGDbkOSe8IYtBZs2bN4OHhgdGjR+P8+fM4c+YMJk6cCC8vL2Vf6ePHj+Hg4IDz58+rvPb27duIjo7GqFGjCu133759+Pfff3Hjxg3cvn0bK1aswM8//4yvv/5aK+9LSMozjUoqlWLt2rX47bff8Msvv8DT0xMymUwbYRJC1IwxhuDgYPTv3x/du3fH2bNnKVlzRaMN7mqUnp7OvL29mampKTMzM2MjR45kmZmZyufv3bvHALDIyEiV1wUFBbF69eqx/Pz8Qvs8dOgQc3R0ZKampszExIS1bt2arVy5sshtS6IP07oqOo3qyJEjzMLCgtnb27OEhATtBk0IqZTXr18zLy8vBoDNnj273J+N+kJbfdgixjgoU6VjMjIyYG5ujqioKHTu3JnrcDQmNFQxCvv9UeKMKfqvSxqKkJSUBE9PT6SmpmLTpk3o06ePtkImhFTQ/fv3MWjQINy6dQuhoaEYPHgw1yHxVnR0NNzc3CCTyTRaNEYQTeJC9+Hc5aQkriOqmIouTNKkSROcO3cOnTp1Qt++fTFjxgzk5uZqI2RCSAXs2rULbdq0QXp6Os6cOUPJmi80ev2uJ0pqEl+7VtFkLJGo3oeEaD9OruXn57OlS5cyAwMD9sknn7B79+5xHRIh5D1v375lEydOZADYF198wV6+fMl1SIJAlc50AC0BqUosFmPq1Kk4ffo0UlNT0aZNG+zatYvrsAghABITE+Hq6orVq1dj+fLl2L59O62rwDOUsDWIloAsmouLCy5fvgx3d3d88cUXmDhxIrKzs7kOixC9tWnTJjg5OeHNmzc4d+4cxo8fT5XLeIgStgbREpDFs7CwQFhYGP7++2/8+++/cHV1RWJiItdhEaJXXr9+DX9/fwwbNgwDBgzAxYsX4ejoyHVYpBiUsDWIloAsmUgkwrhx4xAbG4s3b96gbdu2WLlyJeRyOdehEaLzzp07BycnJ2zbtg0hISFYv349qlWrxnVYpASUsDXI37/kK+wiFg7TS61bt0ZcXBx8fHwwbtw49OjRA/f1ufmBEA3Kzs7GjBkz0LFjR5iZmeHixYvw8/OjJnABoIStQbQEZNmZmppi1apVOHr0KJKSktCyZUusXLmy1KVOCSFld+7cObRp0wZ//PEHFixYgLNnz6JZs2Zch0XKiBK2hlV07rK+6tGjB27cuIGhQ4di3LhxcHd3p6ttohO4rMfw/lV1tWrVcOnSJQQGBsLAQBArLJP/R5XO1EBfKp1p27FjxxAQEICXL19iyZIlGDt2LDXbEUEqai35slQJVIdz585h5MiRuHv3LubOnYtp06bpTaJOSlLM1rl/XzFmyN9f0fKpblTpjOi9gqttb29v5dX2f//9x3VYhJQLV/UYMjMzMXXqVL29qg4JARwcgCVLgLAwxb2Dg6LEslBRwia8ZmZmhn/++QdHjhzB/fv30apVKwQGBiIrK4vr0AgpE23XY2CMYfPmzbC3t1euQHj27Fl8/PHH6j0Qj+lq0SpK2EQQevbsiZs3b+KHH37AH3/8gWbNmiEsLIwGpRHe02Y9hhs3bqBr167w8fFBhw4d8N9//2HmzJl6c1VdQFeLVlHCJoJhbGyM2bNnIz4+Hk5OTvjyyy/Ro0cP3Lp1i+vQCCmWNuoxyGQyfPvtt3B0dERKSgqOHDmCHTt2wM7OrvI7FyBdLVpFCZsITsOGDbFnzx4cOHBA2Uw+ffp0ZGZmch0aIYVosh4DYwwbN26Evb09/vnnH/z000+4du0aevbsWfGd6gBdLVpFCZsIVp8+fXDjxg3MmTMHy5cvh729PVavXo28vDyuQyNESVP1GE6fPo1OnTrB19cXnTt3xn///YfAwEBIpVL1vgEB0tWiVZSwiaAZGxvjhx9+QHx8PNzc3DBmzBi0aNECO3fupP5twhvqrMdw/fp19OvXD59++inevHmDY8eOISwsDPXq1VN32IKlq0WrKGETndCgQQNs2bIFly5dQoMGDTBo0CC4uLjgxIkTXIdGCABFkggOBrZsUdyXN2ncv38fw4cPR+vWrXHr1i1s2bIFcXFxcHd310zAAqeLRav0a+gg0Xlt2rTB4cOHERkZicDAQHTv3h09evRAcHAwnJycuA6vUrRVBEKodPX8PHv2DAsWLMCKFStgaWmJ5cuXIyAgAEZGRlyHxnsFX5J0BiOVJpPJGAAWFRXFdSjkPXK5nO3atYs5ODgwAGzIkCHs5s2bXIdVIWvXMiYWMyaRqN6HhHAdGT/o4vl58eIFmzNnDjM1NWVmZmZswYIFLCsri+uwSBGioqIYACaTyTR6HGoSJzpLJBLh888/x/Xr17FmzRpl8YjPP/8c58+f5zq8MtPVIhDqomvnJyUlBTNmzICdnR0WLVqEsWPH4u7du/juu+9gYmLCdXiEQ5Swic4zMDCAv78/7ty5gzVr1iA+Ph4uLi7o3r07jh8/zvvBabpaBEJddOX83LlzB2PHjkWDBg2watUqTJw4Effv38fSpUtRo0YNrsMjPEAJm+gNIyMj+Pv7Iz4+Htu3b8erV6/Qo0cPODs7Y9euXZDL5VyHWCRdLQKhLkI/P1evXoW3tzeaNm2KPXv2YO7cuXj48CGCg4NhbW3NdXiERyhhE70jkUgwaNAgXLx4EUeOHIGpqSm++OILfPzxx1i7di3evn3LdYgqdLUIhLoI8fwwxhAZGYm+ffvC0dER586dw59//on79+8jMDAQ5ubmXIdIeIgSNtFbIpEIPXv2RGRkJGJiYtC0aVOMGjUKderUwbRp03CbJ52fuloEQl2EdH5kMhn++usvfPzxx+jWrRsePHiADRs2IDExEePHj0eVKlW4DpHwGCVsQgB88skn2Lt3L5KSkjBq1CiEhISgSZMm8PDwQHh4OPLz8zmLTVeLQKiLEM7P1atX8dVXX6FOnTr49ttv0aJFC0RGRuL69esYNmwYDA0NuQ6RCICI8X3EjQBkZGTA3NwcUVFR6Ny5M9fhEDV4+/Yttm/fjr///huxsbGws7PD2LFjERAQwFm/4u3bigRUMM84IIAfyYgv1Hl+1DGnOycnBzt37sTff/+NM2fOwNbWFmPHjsWoUaNga2tbscAIL0VHR8PNzQ0ymQxmZmYaOw4lbDWghK3b4uLisGLFCmzevBl5eXnw9PTEsGHD0Lt3b6rbrINCQhTTxEQiRZN6wf2aNaVXyWKM4eLFi9i0aRM2b96MZ8+eoXv37hg/fjz69etHV9I6SlsJm5rECSmFk5MT/v33Xzx+/BiLFy/GvXv38Pnnn8PGxgZjx45FdHQ0b0eYk/Kp6JzuO3fuYN68eXBwcICzszO2bdsGX19f3Lp1C8ePH8fAgQMpWZNKo4RNSBlZWlpi8uTJiIuLw82bNzF+/HgcOXIEbm5uaNiwIb777jvcvHmT6zAFLykJCAoCvL0V90lJ2jt2eeZ0P3v2DMuXL0eHDh3QuHFjLFmyBK6urjh27BgePXqEX375BQ4ODtoJnOgFahJXA2oS119yuRxnz57Fxo0bERYWhpcvX6J169YYNGgQPD090bJlS4iKywCkkMo0R6uDtzcQFqa4ov6QWAx4er5F376bsHv3bhw9ehQA4OHhAR8fH3h6eqJq1aqaD5LwDvVhCwglbAIoBhkdPnwYmzdvxsGDB5GVlYX69evD09MTnp6e6Ny5My3YUIKkJMDBofhkmZCg+UF2QUHAkiWKZvDC8gAsgVj8Azp06AAvLy8MGTIENWvW1GxQhPeoD5sQgZFKpejfvz+2bduG58+f48iRI/jss8+wZ88e9OjRAzVr1sSXX36JTZs24cWLF1yHyzt8KDGqmNPNAHx4HcMAiLF4sQPS0tJw6tQpTJgwQTDJmstuBqI+tLwmIRoglUrRs2dP9OzZE3/++SeuXr2K8PBwhIeHY9iwYZBIJGjXrh26du2KLl26oGPHjjA1NeU6bE5xVWI0Pz8fly9fxsmTJ3Hy5EkYGdVFdvZyAAwikfj/v0SIsGaNGH5+n2smCA0qqpth8WLtdTMQ9aGETYiGiUQiODo6wtHREbNnz8bjx49x4MABREZGIiQkBAsXLoSBgQHat2+PLl26KBO4vq3MpK0So/n5+bh69SoiIyNx8uRJREdHIyMjA1WqVEHHjh3xww/10KLFf4iJaY4HD0SCnvP+/qj3DwUEAJ06CfN96Svqw1YD6sMmFcUYQ0JCAk6ePKlMIE+fPoWBgQGcnZ3h6uoKJycnODk5oXHjxhCLdbcXS1N92E+fPkVcXBwuXryI8+fP49SpU5DJZDA2NkbHjh3RpUsXdO3aFe3bt9e5MQYl9clLJMD06UBwsPbj0jXa6sOmK2xCOCQSieDg4AAHBwd89dVXYIzhv//+UzbP7ty5E7/88gsAwMzMDG3btlUmcF1L4gUlRgMCih4lXpZk/X5yjouLQ1xcHB49egQAsLCwQLt27TB16lR06dIFzs7OOl/4RugrmRFVgknYCxYswIEDB3DlyhUYGRnh1atXpb6GMYY5c+Zg9erVePXqFTp27IgVK1agyXs1Bl+8eIGvv/4a+/btg1gsxhdffIE//vhD7/sTCTdEIhGaNWuGZs2aYdy4cQCA9PR0ZfKJi4srlMRbtmwJe3t7NG3aFE2bNoW9vT0++ugjQSYjPz9FM21JJUblcjmSk5ORmJiovCUkJODmzZvK5GxpaQknJyf4+Pgov9w0bNhQ76bYCXElM1I8wTSJz5kzBxYWFnj06BHWrFlTpoS9aNEiBAcHY926dWjYsCFmzZqF69evIz4+HsbGxgCA3r17IyUlBatWrUJubi5GjhyJ9u3bY/PmzWWOjZrEibalp6fj0qVLuHjxIm7evKlMWhkZGQAAsViM+vXrK5N4kyZNUKdOHdSuXRu2trawsbHhbUKXy+V49uwZnjx5gpSUFDx58gR3795VJuekpCRkZ2cDAAwNDdG4cWM0bdoUDg4Oep2ci8KHqXL6gOZhFyM0NBSTJ08uNWEzxmBra4upU6di2rRpABRL21lbWyM0NBReXl64desWmjdvjgsXLqBdu3YAgMOHD6NPnz549OhRmQv0U8ImfMAYw9OnT1WuOgv+fefOHbx7905l+xo1aigTeMF99erVYWpqimrVqqFatWpF/tvExAQSiaTYhMgYg1wuR3Z2NjIzM5GZmYmsrKxC/87KyoJMJkNKSooyMaekpCA1NVVldTSRSIS6deuqtCIUtCTY2dnBwEAwDYWcCA0tvpuBRomrB/VhV9K9e/eQmpoKd3d35WPm5uZwcXFBTEwMvLy8EBMTo+zXKuDu7g6xWIzY2Fh8/nnRUzhycnKQk5Oj/FkmkwEArly5opk3Q0g5NWnSRKXrhzEGmUyG9PT0Im8XL15Eeno6MjMz8ebNG5T1e7xIJFLeChJ1WRkbG6Nq1aqoUaMGatSoAWtra3z88ceoUaMGrKyslI9bWloWmZQfPXqkbAInxWvUCNi4UYr9+6sjLc0Q1ta5+OyzF6hTJwfR0VxHpxsKPvs1ff2rswk7NTUVAAothWhtba18LjU1FbVq1VJ53sDAANWrV1duU5Tg4GDMnTu30OOTJk2qbNiECApjrMIfUtnZ2cjOzsaLFy+QRJU8tKocPX6kHNLT02Fubq6x/XOasAMDA7Fo0aISt7l16xbvCugHBQVhypQpyp9fvXqF+vXr4+HDhxr9z9KUjIwM1KtXD8nJyRptztEUip9bFD+3hB4/IPz3IJPJYGdnh+rVq2v0OJwm7KlTp8KvlE6URo0aVWjfNjY2AIC0tDTUrl1b+XhaWhocHR2V2zx9+lTldXl5eXjx4oXy9UWRSqVFDtgxNzcX5C9bATMzM4qfQxQ/tyh+7gn9PWh6iiWnCbtmzZoaq8XbsGFD2NjYICIiQpmgMzIyEBsbq5wu4+rqilevXiEuLg5OTk4AgBMnTkAul8PFxUUjcRFCCCEVIZiKCw8fPsSVK1fw8OFD5Ofn48qVK7hy5QqysrKU2zg4OGD37t0AFINhJk+ejJ9++gnh4eG4fv06hg8fDltbWwwYMAAA0KxZM3h4eGD06NE4f/48zpw5g4kTJ8LLy6vMI8QJIYQQbRDMoLPZs2dj3bp1yp/btGkDAIiMjESXLl0AAAkJCcoR2wAwY8YMvH79GmPGjMGrV6/QqVMnHD58WDkHGwA2bdqEiRMnonv37srCKcuWLStXbFKpFHPmzOHtvNbSUPzcovi5RfFzT+jvQVvxC24eNiGEEKKPBNMkTgghhOgzStiEEEKIAFDCJoQQQgSAEjYhhBAiAJSwy2DBggXo0KEDqlatCgsLizK9hjGG2bNno3bt2qhSpQrc3d0LlV988eIFfHx8YGZmBgsLCwQEBKhMU1OX8h7n/v37KjWi379t375duV1Rz2/dupXz+AGgS5cuhWL76quvVLZ5+PAh+vbti6pVq6JWrVqYPn068vLyOI+/YMlXe3t7VKlSBXZ2dvjmm29UZkAAmj3/y5cvR4MGDWBsbAwXFxecP3++xO23b98OBwcHGBsbo2XLljh48KDK82X5e1Cn8sS/evVqfPrpp7C0tISlpSXc3d0Lbe/n51foXHt4ePAi/tDQ0EKxvT8TBuD3+S/qb1UkEqFv377KbbR5/qOjo9GvXz/Y2tpCJBJhz549pb7m5MmTaNu2LaRSKRo3bozQ0NBC25T3b6pIjJRq9uzZ7Ndff2VTpkxh5ubmZXrNwoULmbm5OduzZw+7evUq8/T0ZA0bNmRv375VbuPh4cFat27Nzp07x06dOsUaN27MvL291R5/eY+Tl5fHUlJSVG5z585lpqamLDMzU7kdABYSEqKy3fvvj6v4GWPMzc2NjR49WiU2mUym8h5btGjB3N3d2eXLl9nBgweZlZUVCwoK4jz+69evs4EDB7Lw8HB2+/ZtFhERwZo0acK++OILle00df63bt3KjIyM2Nq1a9nNmzfZ6NGjmYWFBUtLSyty+zNnzjCJRMIWL17M4uPj2Q8//MAMDQ3Z9evXlduU5e9BXcob/9ChQ9ny5cvZ5cuX2a1bt5ifnx8zNzdnjx49Um4zYsQI5uHhoXKuX7x4ofbYKxJ/SEgIMzMzU4ktNTVVZRs+n//09HSV2G/cuMEkEgkLCQlRbqPN83/w4EH2/fffs127djEAbPfu3SVuf/fuXVa1alU2ZcoUFh8fz/78808mkUjY4cOHlduU95wUhxJ2OYSEhJQpYcvlcmZjY8OWLFmifOzVq1dMKpWyLVu2MMYYi4+PZwDYhQsXlNscOnSIiUQi9vjxY7XFrK7jODo6Mn9/f5XHyvLLXFkVjd/NzY1NmjSp2OcPHjzIxGKxygfbihUrmJmZGcvJyVFL7Iyp7/yHhYUxIyMjlpubq3xMU+ff2dmZTZgwQflzfn4+s7W1ZcHBwUVuP2TIENa3b1+Vx1xcXNjYsWMZY2X7e+Ay/g/l5eWxatWqsXXr1ikfGzFiBOvfv7+6Qy1SeeMv7XNJaOf/t99+Y9WqVWNZWVnKx7R5/t9Xlr+xGTNmsI8//ljlsS+//JL16tVL+XNlz0kBahLXgNKW9gRQ6tKe6qKO48TFxeHKlSsICAgo9NyECRNgZWUFZ2dnrF27Vu3Ly1Um/k2bNsHKygotWrRAUFAQ3rx5o7Lfli1bqqzm1qtXL2RkZODmzZu8iP99BevsfrjMpLrP/7t37xAXF6fyuysWi+Hu7q783f1QTEyMyvaA4lwWbF+Wvwd1qUj8H3rz5g1yc3MLLeRw8uRJ1KpVC/b29hg3bhzS09PVGjtQ8fizsrJQv3591KtXD/3791f5HRba+V+zZg28vLxgYmKi8rg2zn9FlPb7r45zUkAwlc6ERJNLe1YklsoeZ82aNWjWrBk6dOig8vi8efPQrVs3VK1aFUePHsX48eORlZWFb775hvP4hw4divr168PW1hbXrl3DzJkzkZCQgF27din3W9T/T8FzXMf/vufPn2P+/PkYM2aMyuOaOP/Pnz9Hfn5+kefmv//+K/I1xZ3L93/XCx4rbht1qUj8H5o5cyZsbW1VPmA9PDwwcOBANGzYEHfu3MF3332H3r17IyYmBhKJhNP47e3tsXbtWrRq1QoymQxLly5Fhw4dcPPmTdStW1dQ5//8+fO4ceMG1qxZo/K4ts5/RRT3+5+RkYG3b9/i5cuXlf6dLKC3CVuoS3sWKGv8lfX27Vts3rwZs2bNKvTc+4+1adMGr1+/xpIlS8qUMDQd//vJrWXLlqhduza6d++OO3fu4KOPPqrwfgto6/xnZGSgb9++aN68OX788UeV5ypz/knRFi5ciK1bt+LkyZMqA7e8vLyU/27ZsiVatWqFjz76CCdPnkT37t25CFXJ1dUVrq6uyp87dOiAZs2aYdWqVZg/fz6HkZXfmjVr0LJlSzg7O6s8zufzr016m7CFurRngbLGX9nj7NixA2/evMHw4cNL3dbFxQXz589HTk5OqTV1tRX/+7EBwO3bt/HRRx/Bxsam0CjNtLQ0AODN+c/MzISHhweqVauG3bt3w9DQsMTty3P+i2NlZQWJRKI8FwXS0tKKjdfGxqbE7cvy96AuFYm/wNKlS7Fw4UIcP34crVq1KnHbRo0awcrKCrdv31ZrwqhM/AUMDQ3Rpk0b3L59G4Bwzv/r16+xdetWzJs3r9TjaOr8V0Rxv/9mZmaoUqUKJBJJpf9PlcrV463nyjvobOnSpcrHZDJZkYPOLl68qNzmyJEjGht0VtHjuLm5FRqdXJyffvqJWVpaVjjWoqjrPJ0+fZoBYFevXmWM/W/Q2fujNFetWsXMzMxYdnY25/HLZDL2ySefMDc3N/b69esyHUtd59/Z2ZlNnDhR+XN+fj6rU6dOiYPOPvvsM5XHXF1dCw06K+nvQZ3KGz9jjC1atIiZmZmxmJiYMh0jOTmZiUQitnfv3krH+6GKxP++vLw8Zm9vz7799lvGmDDOP2OKz1epVMqeP39e6jE0ef7fhzIOOmvRooXKY97e3oUGnVXm/1QZT7m21lMPHjxgly9fVk5tunz5Mrt8+bLKFCd7e3u2a9cu5c8LFy5kFhYWbO/evezatWusf//+RU7ratOmDYuNjWWnT59mTZo00di0rpKO8+jRI2Zvb89iY2NVXpeUlMREIhE7dOhQoX2Gh4ez1atXs+vXr7OkpCT2999/s6pVq7LZs2dzHv/t27fZvHnz2MWLF9m9e/fY3r17WaNGjVjnzp2VrymY1tWzZ0925coVdvjwYVazZk2NTesqT/wymYy5uLiwli1bstu3b6tMZcnLy2OMafb8b926lUmlUhYaGsri4+PZmDFjmIWFhXJEva+vLwsMDFRuf+bMGWZgYMCWLl3Kbt26xebMmVPktK7S/h7UpbzxL1y4kBkZGbEdO3aonOuCv+/MzEw2bdo0FhMTw+7du8eOHz/O2rZty5o0aaLWL3cVjX/u3LnsyJEj7M6dOywuLo55eXkxY2NjdvPmTZX3yNfzX6BTp07syy+/LPS4ts9/Zmam8jMeAPv111/Z5cuX2YMHDxhjjAUGBjJfX1/l9gXTuqZPn85u3brFli9fXuS0rpLOSVlRwi6DESNGMACFbpGRkcpt8P9zYgvI5XI2a9YsZm1tzaRSKevevTtLSEhQ2W96ejrz9vZmpqamzMzMjI0cOVLlS4C6lHace/fuFXo/jDEWFBTE6tWrx/Lz8wvt89ChQ8zR0ZGZmpoyExMT1rp1a7Zy5coit9V2/A8fPmSdO3dm1atXZ1KplDVu3JhNnz5dZR42Y4zdv3+f9e7dm1WpUoVZWVmxqVOnqkyb4ir+yMjIIn/fALB79+4xxjR//v/8809mZ2fHjIyMmLOzMzt37pzyOTc3NzZixAiV7cPCwljTpk2ZkZER+/jjj9mBAwdUni/L34M6lSf++vXrF3mu58yZwxhj7M2bN6xnz56sZs2azNDQkNWvX5+NHj263B+2mop/8uTJym2tra1Znz592KVLl1T2x+fzzxhj//33HwPAjh49Wmhf2j7/xf39FcQ8YsQI5ubmVug1jo6OzMjIiDVq1EglFxQo6ZyUFS2vSQghhAgAzcMmhBBCBIASNiGEECIAlLAJIYQQAaCETQghhAgAJWxCCCFEAChhE0IIIQJACZsQQggRAErYhBBCiABQwiaEEEIEgBI2IYQQIgCUsAkhhBABoIRNCCm3Z8+ewcbGBj///LPysbNnz8LIyAgREREcRkaI7qLFPwghFXLw4EEMGDAAZ8+ehb29PRwdHdG/f3/8+uuvXIdGiE6ihE0IqbAJEybg+PHjaNeuHa5fv44LFy5AKpVyHRYhOokSNiGkwt6+fYsWLVogOTkZcXFxaNmyJdchEaKzqA+bEFJhd+7cwZMnTyCXy3H//n2uwyFEp9EVNiGkQt69ewdnZ2c4OjrC3t4ev//+O65fv45atWpxHRohOokSNiGkQqZPn44dO3bg6tWrMDU1hZubG8zNzbF//36uQyNEJ1GTOCGk3E6ePInff/8dGzZsgJmZGcRiMTZs2IBTp05hxYoVXIdHiE6iK2xCCCFEAOgKmxBCCBEAStiEEEKIAFDCJoQQQgSAEjYhhBAiAJSwCSGEEAGghE0IIYQIACVsQgghRAAoYRNCCCECQAmbEEIIEQBK2IQQQogAUMImhBBCBOD/AEV0qB0BEDeDAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.set_xlim(-1, 1)\n",
    "ax.set_ylim(-1, 1)\n",
    "ax.add_patch(plt.Rectangle((-1,-1), 2, 2, fill=False))\n",
    "ax.add_artist(plt.Circle((0, 0), 1, fill=False))\n",
    "x = np.random.uniform(-1, 1, size=100)\n",
    "y = np.random.uniform(-1, 1, size=100)\n",
    "ax.scatter(x, y, s=25, c='blue', marker='o')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "By laws of large numbers, $\\pi$ can be approximated by a stochastic algorithm\n",
    "\n",
    "$$\n",
    "E\\left[\\boldsymbol{1}\\left\\{ x^{2}+y^{2}\\leq1\\right\\} \\right] = \\frac{\\pi r^{2}}{\\left(2r\\right)^{2}}= \\frac{\\pi}{4}\n",
    "$$\n",
    "\n",
    "it implies  $\\pi=4\\times E[ 1 \\{  x^{2}+y^{2}\\leq1 \\}]$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The estimated pi = 3.142011\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "n = 20000000\n",
    "Z = 2 * np.random.rand(n, 2) - 1\n",
    "\n",
    "inside = np.mean(np.sqrt(np.sum(Z**2, axis=1)) <= 1)\n",
    "print(f\"The estimated pi = {inside * 4:.6f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Sample size can be made as large as the computer's memory permits.\n",
    "* Iterate it with average of averages and so on, for higher accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Writing Script\n",
    "\n",
    "* A script is a piece of code for a particular purpose. Thousands of lines are not written from the beginning to the end. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Recursively development\n",
    "  * small manageable tasks\n",
    "  * test code constantly \n",
    "  * encapsulate into DIY functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Integrate small pieces into the super structure. \n",
    "* Add comments to the script to facilitate readability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Finite Sample Evaluation\n",
    "\n",
    "\n",
    "* Real world sample is finite\n",
    "* Asymptotic theory is a mathematical apparatus to approximate finite sample distributions\n",
    "* Modern econometric theory is built on asymptotics\n",
    "* Simulation is one way to evaluate the accuracy of approximation.\n",
    "\n",
    "## Example\n",
    "\n",
    "For the OLS estimator, \n",
    "\n",
    "* Classical views $X$ as fixed regressions and only cares about the randomness of the error term.\n",
    "* Modern econometrics textbook emphasizes that a random $X$ is more appropriate\n",
    "for econometrics applications. \n",
    "* In rigorous textbooks, the moment of $X$ is explicitly stated as $E[X_i X_i'] < \\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* A Pareto distribution with shape coefficient between 1 and 2 has finite population mean, but infinite variance. \n",
    "* If $X$ follows a\n",
    "[Pareto distribution](https://en.wikipedia.org/wiki/Pareto_distribution) with shape coefficient 1.5, it violates the assumptions for OLS stated in most of econometric textbooks.\n",
    "* Question: Is asymptotic inferential theory for the OLS estimator valid? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We write a script to investigate this problem. The following steps develop the code.\n",
    "\n",
    " 1. given a sample size, get the OLS `b_hat` and its associated `t_value`.\n",
    " 2. wrap `t_value` as a user-defined function so that we can reuse it for many times.\n",
    " 3. given a sample size, report the size under two distributions.\n",
    " 4. wrap step 3 again as a user-defined function, ready for different sample sizes.\n",
    " 5. develop the super structure to connect the workhorse functions.\n",
    " 6. add comments and documentation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### $t$-statistic\n",
    "\n",
    "$$\n",
    "\\sqrt{n} (\\hat{\\beta} - \\beta_0) |X = (X'X/n)^{-1}  (X' e /\\sqrt{n}),\n",
    "$$\n",
    "\n",
    "the $k$-th element of the vector coefficient conditional on $X$ is\n",
    "$$\n",
    "\\widehat{\\beta}_{k}|X=\\eta_{k}'\\widehat{\\beta}|X\n",
    "\\sim N\\left(\\beta_{k},\\sigma^{2}\\left(X'X\\right)_{kk}^{-1}\\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm, t\n",
    "\n",
    "def simulation(n, dist='normal', df=5):\n",
    "    if dist == 'normal':\n",
    "        e = np.random.normal(size=n)\n",
    "    elif dist == 't':\n",
    "        e = np.random.standard_t(df, size=n)\n",
    "\n",
    "    b0 = np.array([1, 2])\n",
    "    X = np.column_stack((np.ones(n), np.random.pareto(1.5, size=n)))\n",
    "    Y = X.dot(b0) + e\n",
    "\n",
    "    bhat = np.linalg.inv(X.T @ X) @ (X.T @ Y)\n",
    "    bhat2 = bhat[1] # parameter we want to test\n",
    "\n",
    "    e_hat = Y - X.dot(bhat)\n",
    "    sigma_hat_square = np.sum(e_hat**2) / (n - 2)\n",
    "    sig_B = np.linalg.inv(X.T.dot(X)) * sigma_hat_square\n",
    "    t_value_2 = (bhat2 - b0[1]) / np.sqrt(sig_B[1, 1])\n",
    "\n",
    "    return t_value_2\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.stats import t, norm\n",
    "\n",
    "def report(n):\n",
    "    # collect the test size from the two distributions\n",
    "    # this function contains some repetitive code, but is OK for such a simple one\n",
    "    TEST_SIZE = np.zeros(3)\n",
    "\n",
    "    # e ~ normal distribution, under which the t-dist is exact\n",
    "    Res = np.array([simulation(n, 'normal') for i in range(Rep)])\n",
    "    TEST_SIZE[0] = np.mean(np.abs(Res) > t.ppf(0.975, n-2))\n",
    "    TEST_SIZE[1] = np.mean(np.abs(Res) > norm.ppf(0.975))\n",
    "\n",
    "    # e ~ t-distribution, under which the exact distribution is complicated.\n",
    "    # we rely on asymptotic normal distribution for inference instead\n",
    "    Res = np.array([simulation(n, 't', df) for i in range(Rep)])\n",
    "    TEST_SIZE[2] = np.mean(np.abs(Res) > norm.ppf(0.975))\n",
    "\n",
    "    return TEST_SIZE\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       n  exact  normal.asym  cauchy.asym\n",
      "0    5.0  0.049        0.144        0.150\n",
      "1   10.0  0.052        0.078        0.102\n",
      "2  200.0  0.050        0.050        0.030\n"
     ]
    }
   ],
   "source": [
    "Rep = 1000\n",
    "df = 1  # t dist. with df = 1 is Cauchy\n",
    "\n",
    "# run the calculation of the empirical sizes for different sample sizes\n",
    "NN = [5, 10, 200]\n",
    "RES = np.zeros((len(NN), 4))\n",
    "RES[:, 0] = NN\n",
    "for i, n in enumerate(NN):\n",
    "    test_sizes = report(n)\n",
    "    RES[i, 1] = test_sizes[0]  # exact\n",
    "    RES[i, 2] = test_sizes[1]  # normal asymptotic\n",
    "    RES[i, 3] = test_sizes[2]  # Cauchy asymptotic\n",
    "\n",
    "print(pd.DataFrame(RES, columns=['n', 'exact', 'normal.asym', 'cauchy.asym']))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Simulation Results\n",
    "\n",
    "* 1st column: the error is normal, use exact distribution theory to find the critical value (according to $t$-distribution.)\n",
    "\n",
    "* 2nd column: still uses the normal error\n",
    "  * change the critical value to asymptotic normal distribution\n",
    "  \n",
    "* 3rd column: error distribution is Cauchy\n",
    "  * asymptotic approximation breaks down\n",
    "  * test size does not converge to the nominal 5%  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Observations\n",
    "\n",
    "* $X$ is always Pareto \n",
    "* distribution of $X$ doesn't matter in our simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Justification\n",
    "\n",
    "* Self-normalized $t$ statistic does not break down despite that $X'X/n$ does not converge\n",
    "* Regardless the distribution of $X$, when the error term is normal \n",
    "  * numerator follows $N(0,1)$\n",
    "  * demonimator follows $\\chi^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### World View\n",
    "\n",
    "\n",
    "* Fundamental question: how to quantify uncertainty.\n",
    "  * data $(X_1,X_2,\\ldots,X_n)$\n",
    "  * sample mean $\\bar{X}$\n",
    "  * sample variance $s$\n",
    "  * frequentist confidence interval about the population mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Asymptotics is imaginary.\n",
    "* Let's be realistic: we have a finite sample $n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bootstrap\n",
    "\n",
    "* Let $X_1, X_2, \\ldots, X_n \\sim F$ be an i.i.d. sample of $n$ observations following a distribution $F$. \n",
    "* The finite sample distribution of a statistic $T_n(\\theta)\\sim G_n(\\cdot, F)$ usually depends on the sample size $n$, as well as the unknown true distribution $F$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Key Idea\n",
    "\n",
    "* Bootstrap replaces the unknown distribution $F$ in $G_n(\\cdot, F)$ by the empirical distribution function\n",
    "\n",
    "$$\n",
    "\\hat{F}_n(\\cdot) = n^{-1} \\sum_{i=1}^n 1\\{\\cdot \\leq X_i\\}\n",
    "$$\n",
    "\n",
    "* Bootstrap inference is drawn from the bootstrap distribution\n",
    "\n",
    "$$\n",
    "G_n(\\cdot, \\hat{F}_n)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Key Idea\n",
    "\n",
    "* Bootstrap replaces the unknown distribution $F$ in $G_n(\\cdot, F)$ by the empirical distribution function\n",
    "\n",
    "$$\n",
    "\\hat{F}_n(\\cdot) = n^{-1} \\sum_{i=1}^n 1\\{\\cdot \\leq X_i\\}\n",
    "$$\n",
    "\n",
    "* Bootstrap inference is drawn from the bootstrap distribution\n",
    "\n",
    "$$\n",
    "G_n(\\cdot, \\hat{F}_n)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Compare to Asymptotic Theory\n",
    "\n",
    "* Bootstrap is a finite-sample practice \n",
    "* It doesn't refer to an imaginary world where $n\\to \\infty$ at its face value\n",
    "\n",
    "\n",
    "* Asymptotic theory approximates $G_n(\\cdot, F)$ by its limit \n",
    "\n",
    "$$G(\\cdot, F) := \\lim_{n\\to\\infty} G_n(\\cdot, F).$$ \n",
    "\n",
    "* In many cases $G(\\cdot, F)$ is independent of $F$ and it becomes $G(\\cdot)$. Such a $T_n(\\theta)$ is called *asymptotically pivotal*, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGzCAYAAAAIWpzfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAt40lEQVR4nO3df1iUdb7/8dcwMgOooIWiEkHZZj9MbSVZLLM2WvqFWbm52SK6mcfSSqlN/BGUZnQqjXNK41LTPP1YLU9ZJ13NSL9tSXki2W0zbc1fnRKUTcXEGBk+3z+6mG0ClEFg+MDzcV1zXc5nPp/7ft/zGZzX3Pc99ziMMUYAAAAWCAl2AQAAAA1FcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAVpQVVWVHnzwQcXFxSkkJETDhw8PdkkBu+KKK3TFFVf4tZWWlmrEiBE6/fTT5XA4lJeXF5TaTsXu3bvlcDj0wgsv+NoefvhhORyOFln/z5/XjRs3yuFwaOXKlS2y/jFjxighIaFF1gWcCoILrPDCCy/I4XDUe/voo4/8+v/www96+umnlZSUpKioKIWFhencc8/VpEmT9OWXX/r61bwx1dwiIiJ05plnKi0tTUuXLlVlZWWtWsaMGVNvHWvXrj3hdixZskRPPvmkRowYoWXLlmnKlClN8wQF2ZQpU7Ru3TpNmzZNL774oq655ppglxQ03377rR5++GEVFxcHu5RaWnNtQEN1CHYBQCBmzZqls846q1b7Oeec4/t3WVmZrrnmGhUVFemGG27QqFGj1KlTJ23fvl3Lly/XwoUL5fF4/MY/99xz6tSpkyorK/XNN99o3bp1+sMf/qC8vDy9/fbbiouL8+vvdru1ePHiWnX079//hPW/9957io2N1dNPPx3IZrcq77zzTq229957TzfeeKMeeOCBIFTUfGbOnKmsrKyAxnz77bd65JFHlJCQoAEDBjR4XF3Pa1M7UW2LFi1SdXV1s9cAnCqCC6xy7bXXKjEx8YR9xowZoy1btmjlypW65ZZb/B6bPXu2ZsyYUWvMiBEjFB0d7bufnZ2tl19+WaNHj9Zvf/vbWnt0OnTooN///vcB179//3516dIl4HH1qa6ulsfjUVhYWJMt82RcLlettqberh9++EEul0shIcHdKdyhQwd16NC8/01WVFQoIiKizue1JYWGhgZ1/UBDcagIbcrHH3+s1atX64477qgVWqQf95Q89dRTDVrW7bffrnHjxunjjz/W+vXrT6mumvMnNmzYoM8//9x3aGnjxo2SpKNHj+r+++9XXFyc3G63+vTpo6eeeko///F2h8OhSZMm6eWXX9aFF14ot9t9wsNTDodDDz/8cK32hIQEjRkzxne/5lDchx9+qMzMTHXr1k0dO3bUTTfdpAMHDviN/em5GDXjjDGaP3++b7tq7Ny5U7/97W912mmnKSIiQr/61a+0evVqv+XVnMuxfPlyzZw5U7GxsYqIiFB5ebnGjBmjTp06ae/evbrhhhvUqVMnxcbGav78+ZKkzz77TL/+9a/VsWNHxcfH65VXXjnZVEiSDh06pDFjxigqKkpdunRRRkaGDh06VKtfXee4rF+/Xpdddpm6dOmiTp06qU+fPpo+fbpvWy655BJJ0tixY33PR815M1dccYX69u2roqIiXX755YqIiPCNrevcIUnyer2aPn26evTooY4dO2rYsGH6+uuv/fr8fD5r/HSZJ6utrnNcAn1drlq1Sn379pXb7daFF1540kOnQGOwxwVWOXz4sMrKyvzaHA6HTj/9dEnSW2+9JUlKT09vkvWlp6dr4cKFeuedd3T11Vf7PfbzOkJDQxUVFVXncrp166YXX3xRc+bM0ffff6/c3FxJ0vnnny9jjIYNG6YNGzbojjvu0IABA7Ru3Tr98Y9/1DfffFPrsNJ7772nV199VZMmTVJ0dHSTnlB5zz33qGvXrsrJydHu3buVl5enSZMmacWKFXX2v/zyy/Xiiy8qPT1dV199tUaPHu17rLS0VIMHD1ZFRYXuvfdenX766Vq2bJmGDRumlStX6qabbvJb1uzZs+VyufTAAw+osrLStwfC6/Xq2muv1eWXX64nnnhCL7/8siZNmqSOHTtqxowZuv3223XzzTcrPz9fo0ePVnJycp2HE2sYY3TjjTfqgw8+0IQJE3T++efrjTfeUEZGxkmfn88//1w33HCD+vXrp1mzZsntdmvHjh368MMPJf04n7NmzVJ2drbGjx+vIUOGSJIGDx7sW8Y///lPXXvttfrd736n3//+94qJiTnhOufMmSOHw6GpU6dq//79ysvLU0pKioqLixUeHn7Smms0pLafCvR1+cEHH+j111/X3Xffrc6dO+s///M/dcstt2jv3r2+v0+gSRjAAkuXLjWS6ry53W5fv5tuuslIMgcPHmzQcnNycowkc+DAgTofP3jwoJFkbrrpJl9bRkZGnXUMHTr0pOsbOnSoufDCC/3aVq1aZSSZRx991K99xIgRxuFwmB07dvjaJJmQkBDz+eefN2j7JJmcnJxa7fHx8SYjI8N3v+b5TUlJMdXV1b72KVOmGKfTaQ4dOuS3DT/fVklm4sSJfm2TJ082ksxf/vIXX9uRI0fMWWedZRISEozX6zXGGLNhwwYjyZx99tmmoqLCbxk1z/Vjjz3mazt48KAJDw83DofDLF++3Ne+bdu2erf3p2qe7yeeeMLXVlVVZYYMGWIkmaVLl/raa14fNZ5++ukTvl6MMeZ///d/ay2nxtChQ40kk5+fX+djP31ea56X2NhYU15e7mt/9dVXjSTzH//xH762n89nfcs8UW0ZGRkmPj7edz/Q16XL5fJr++tf/2okmWeeeabWuoBTwaEiWGX+/Plav3693+3Pf/6z7/Hy8nJJUufOnZtkfZ06dZIkHTlyxK89LCysVh1z585t1DrWrFkjp9Ope++916/9/vvvlzHGb/skaejQobrgggsata6TGT9+vN+hkSFDhsjr9WrPnj0BL2vNmjUaNGiQLrvsMl9bp06dNH78eO3evVtbt27165+RkVHvHoRx48b5/t2lSxf16dNHHTt21K233upr79Onj7p06aKdO3eetK4OHTrorrvu8rU5nU7dc889J92mmvN43nzzzUafyOp2uzV27NgG9x89erTf63nEiBHq2bOn1qxZ06j1N1Sgr8uUlBT17t3bd79fv36KjIw86XwAgeJQEawyaNCgE56cGxkZKenHoNEUJ4t+//33kmoHIafTqZSUlFNeviTt2bNHvXr1qrWO888/3/f4T53oMMipOvPMM/3ud+3aVZJ08ODBgJe1Z88eJSUl1Wr/6Xb17dvX117fdoWFhalbt25+bVFRUTrjjDNqnX8SFRV10lr37Nmjnj17+kJpjT59+pxwnCSNHDlSixcv1rhx45SVlaWrrrpKN998s0aMGNHgE4ljY2MDOhH3F7/4hd99h8Ohc845R7t3727wMhoj0Nflz1870o+vn8a8doATYY8L2pTzzjtP0o8nbTaFv//975L8v24dbIGc11Afr9dbZ7vT6ayz3fzsZMzmUN921VdTMGoNDw/X+++/r3fffVfp6en629/+ppEjR+rqq6+u9zmtaxlNrb6L5DW0pqYQzNcO2heCC9qUtLQ0SdJLL73UJMt78cUXJUmpqalNsry6xMfH69tvv611OGrbtm2+xxura9eutb4t4/F4tG/fvkYvs6Hi4+O1ffv2Wu1NsV2nIj4+Xvv27fPtTatRV611CQkJ0VVXXaV58+Zp69atmjNnjt577z1t2LBBUv0horH+8Y9/+N03xmjHjh1+J2XXNc9S7b0igdTWnK9L4FQQXNCmJCcn65prrtHixYu1atWqWo97PJ4GXyTtlVde0eLFi5WcnKyrrrqqiSv9l+uuu05er1fPPvusX/vTTz8th8Oha6+9ttHL7t27t95//32/toULF7bIJ/HrrrtOmzdvVmFhoa/t6NGjWrhwoRISEprtPJ2G1FVVVaXnnnvO1+b1evXMM8+cdOx3331Xq63mQm41V1nu2LGjJNUZJBrjv/7rv/zCw8qVK7Vv3z6/10Xv3r310Ucf+V1Y8e233671telAamvO1yVwKjjHBVb585//7PvE91ODBw/W2WefLenH/+h/85vf6Oabb1ZaWpquuuoqdezYUf/4xz+0fPly7du3r9a1XFauXKlOnTrJ4/H4rpz74Ycfqn///nrttdeadZvS0tJ05ZVXasaMGdq9e7f69++vd955R2+++aYmT57sd8JjoMaNG6cJEybolltu0dVXX62//vWvWrdund/F9ppLVlaW/vSnP+naa6/Vvffeq9NOO03Lli3Trl279N///d9Bu7hcWlqaLr30UmVlZWn37t264IIL9Prrr+vw4cMnHTtr1iy9//77uv766xUfH6/9+/drwYIFOuOMM3wnIffu3VtdunRRfn6+OnfurI4dOyopKanR5yaddtppuuyyyzR27FiVlpYqLy9P55xzju68805fn3HjxmnlypW65pprdOutt+qrr77SSy+9VOu1E0htzfm6BE4FwQVWyc7OrrN96dKlvuDSrVs3bdq0SQsWLNCKFSs0Y8YMeTwexcfHa9iwYbrvvvtqja/5hklYWJiio6M1YMAALVmyRKNGjZLb7W6+DdKPhx7eeustZWdna8WKFVq6dKkSEhL05JNP6v777z+lZd95553atWuXnn/+ea1du1ZDhgzR+vXrm3UPUo2YmBht2rRJU6dO1TPPPKMffvhB/fr10//8z//o+uuvb/b116fm+Z48ebJeeuklORwODRs2THPnztXFF198wrHDhg3T7t27tWTJEpWVlSk6OlpDhw7VI4884ruGT2hoqJYtW6Zp06ZpwoQJqqqq0tKlSxsdXKZPn66//e1vys3N1ZEjR3TVVVdpwYIFioiI8PVJTU3V3LlzNW/ePE2ePFmJiYl6++23a71+AqmtOV+XwKlwGM6cAgAAluAcFwAAYA2CCwAAsAbBBQAAWCPg4PL+++8rLS1NvXr1ksPhqPMrpz+3ceNG/fKXv5Tb7dY555zj+zVSAACAQAQcXI4ePar+/fv7flb+ZHbt2qXrr79eV155pYqLizV58mSNGzdO69atC7hYAADQvp3St4ocDofeeOMNDR8+vN4+U6dO1erVq32XTpek3/3udzp06JDWrl3b2FUDAIB2qNmv41JYWFjrx+hSU1M1efLkesdUVlb6rkIpSdXV1fruu+90+umnN/nltAEAQPMwxujIkSPq1atXk110stmDS0lJiWJiYvzaYmJiVF5ermPHjtX5g2O5ubl65JFHmrs0AADQAr7++mudccYZTbKsVnnl3GnTpikzM9N3//DhwzrzzDP19ddfKzIyMoiVAQBaQoWnSoPmFEiS/t8fr1C4q+5fn0bwHPNUKSzUWeeRkBCHQ2GhTpWXlysuLk6dO3dusvU2e3Dp0aOHSktL/dpKS0sVGRlZ78+7u93uOi+zHhkZSXABgHagg6dKIe4ff9YgJrqrIlyt8nM2GqgpT/No9uu4JCcnq6CgwK9t/fr1Sk5Obu5VAwCANibg4PL999+ruLhYxcXFkn78unNxcbH27t0r6cfDPKNHj/b1nzBhgnbu3KkHH3xQ27Zt04IFC/Tqq69qypQpTbMFAACg3Qg4uHzyySe6+OKLfb+impmZqYsvvtj3q7379u3zhRhJOuuss7R69WqtX79e/fv319y5c7V48WKlpqY20SYAAID2IuCDhldccYVOdOmXuq6Ke8UVV2jLli2BrgoAAMAPv1UEAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1OgS7AABA+2KM0bHj3hP2qfCc+HG0XwQXAECLMcZoRH6hivYcDHYpsBSHigAALebYcW9AoSUxvqvCQ53NWBFswx4XAEBQfDIzRRGuE4eS8FCnHA5HC1UEGxBcAABBEeFyKsLF2xACw6EiAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrNCq4zJ8/XwkJCQoLC1NSUpI2b958wv55eXnq06ePwsPDFRcXpylTpuiHH35oVMEAAKD9Cji4rFixQpmZmcrJydGnn36q/v37KzU1Vfv376+z/yuvvKKsrCzl5OToiy++0PPPP68VK1Zo+vTpp1w8AABoXzoEOmDevHm68847NXbsWElSfn6+Vq9erSVLligrK6tW/02bNunSSy/VqFGjJEkJCQm67bbb9PHHH59i6QAQOGOMjh33BruMdqvCw3OPUxNQcPF4PCoqKtK0adN8bSEhIUpJSVFhYWGdYwYPHqyXXnpJmzdv1qBBg7Rz506tWbNG6enp9a6nsrJSlZWVvvvl5eWBlAkAdTLGaER+oYr2HAx2KQAaKaDgUlZWJq/Xq5iYGL/2mJgYbdu2rc4xo0aNUllZmS677DIZY1RVVaUJEyac8FBRbm6uHnnkkUBKA4CTOnbcS2hpJRLjuyo81BnsMmChgA8VBWrjxo167LHHtGDBAiUlJWnHjh267777NHv2bD300EN1jpk2bZoyMzN998vLyxUXF9fcpQJoRz6ZmaIIF2+cwRIe6pTD4Qh2GbBQQMElOjpaTqdTpaWlfu2lpaXq0aNHnWMeeughpaena9y4cZKkiy66SEePHtX48eM1Y8YMhYTUPj/Y7XbL7XYHUhoABCTC5VSEq9k/uwFoYgF9q8jlcmngwIEqKCjwtVVXV6ugoEDJycl1jqmoqKgVTpzOHz/lGGMCrRcAALRjAX/cyMzMVEZGhhITEzVo0CDl5eXp6NGjvm8ZjR49WrGxscrNzZUkpaWlad68ebr44ot9h4oeeughpaWl+QIMAABAQwQcXEaOHKkDBw4oOztbJSUlGjBggNauXes7YXfv3r1+e1hmzpwph8OhmTNn6ptvvlG3bt2UlpamOXPmNN1WAACAdsFhLDheU15erqioKB0+fFiRkZHBLgeApSo8Vboge50kaeusVM5xAZpZc7x/81tFAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYo0OwCwDQNhljdOy4N9hl+KnwtK56AASO4AKgyRljNCK/UEV7Dga7FABtDIeKADS5Y8e9rTq0JMZ3VXioM9hlAGgE9rgAaFafzExRhKt1hYTwUKccDkewywDQCAQXAM0qwuVUhIv/agA0DQ4VAQAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYo1HBZf78+UpISFBYWJiSkpK0efPmE/Y/dOiQJk6cqJ49e8rtduvcc8/VmjVrGlUwAABovzoEOmDFihXKzMxUfn6+kpKSlJeXp9TUVG3fvl3du3ev1d/j8ejqq69W9+7dtXLlSsXGxmrPnj3q0qVLU9QPAADakYCDy7x583TnnXdq7NixkqT8/HytXr1aS5YsUVZWVq3+S5Ys0XfffadNmzYpNDRUkpSQkHBqVQNoFsYYHTvuPeXlVHhOfRkAUJeAgovH41FRUZGmTZvmawsJCVFKSooKCwvrHPPWW28pOTlZEydO1Jtvvqlu3bpp1KhRmjp1qpxOZ51jKisrVVlZ6btfXl4eSJkAGsEYoxH5hSraczDYpQBAvQI6x6WsrExer1cxMTF+7TExMSopKalzzM6dO7Vy5Up5vV6tWbNGDz30kObOnatHH3203vXk5uYqKirKd4uLiwukTACNcOy4t8lDS2J8V4WH1v0BBQAaI+BDRYGqrq5W9+7dtXDhQjmdTg0cOFDffPONnnzySeXk5NQ5Ztq0acrMzPTdLy8vJ7wALeiTmSmKcJ164AgPdcrhcDRBRQDwo4CCS3R0tJxOp0pLS/3aS0tL1aNHjzrH9OzZU6GhoX6Hhc4//3yVlJTI4/HI5XLVGuN2u+V2uwMpDUATinA5FeFq9s81ABCwgA4VuVwuDRw4UAUFBb626upqFRQUKDk5uc4xl156qXbs2KHq6mpf25dffqmePXvWGVoAAADqE/B1XDIzM7Vo0SItW7ZMX3zxhe666y4dPXrU9y2j0aNH+528e9ddd+m7777Tfffdpy+//FKrV6/WY489pokTJzbdVgAAgHYh4H3BI0eO1IEDB5Sdna2SkhINGDBAa9eu9Z2wu3fvXoWE/CsPxcXFad26dZoyZYr69eun2NhY3XfffZo6dWrTbQUAAGgXHMYYE+wiTqa8vFxRUVE6fPiwIiMjg10O0CZVeKp0QfY6SdLWWamc4wLglDXH+ze/VQQAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGt0CHYBAGozxujYcW+LrrPC07LrA4DGILgArYwxRiPyC1W052CwSwGAVodDRUArc+y4N6ihJTG+q8JDnUFbPwCcCHtcgFbsk5kpinC1bIgID3XK4XC06DoBoKEILkArFuFyKsLFnykA1OBQEQAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGo0KLvPnz1dCQoLCwsKUlJSkzZs3N2jc8uXL5XA4NHz48MasFgAAtHMdAh2wYsUKZWZmKj8/X0lJScrLy1Nqaqq2b9+u7t271ztu9+7deuCBBzRkyJBTKhgIBmOMjh33tsi6Kjwtsx4AsJHDGGMCGZCUlKRLLrlEzz77rCSpurpacXFxuueee5SVlVXnGK/Xq8svv1x/+MMf9Je//EWHDh3SqlWr6l1HZWWlKisrfffLy8sVFxenw4cPKzIyMpBygVNmjNGI/EIV7TnY4uveOitVEa6AP18AQKtQXl6uqKioJn3/DuhQkcfjUVFRkVJSUv61gJAQpaSkqLCwsN5xs2bNUvfu3XXHHXc0aD25ubmKiory3eLi4gIpE2hSx457gxJaEuO7KjzU2eLrBYDWLKCPcmVlZfJ6vYqJifFrj4mJ0bZt2+oc88EHH+j5559XcXFxg9czbdo0ZWZm+u7X7HEBgu2TmSmKcLVMmAgPdcrhcLTIugDAFs26D/rIkSNKT0/XokWLFB0d3eBxbrdbbre7GSsDGifC5eTQDQAEUUD/A0dHR8vpdKq0tNSvvbS0VD169KjV/6uvvtLu3buVlpbma6uurv5xxR06aPv27erdu3dj6gYAAO1QQOe4uFwuDRw4UAUFBb626upqFRQUKDk5uVb/8847T5999pmKi4t9t2HDhunKK69UcXExh38AAEBAAt7nnZmZqYyMDCUmJmrQoEHKy8vT0aNHNXbsWEnS6NGjFRsbq9zcXIWFhalv375+47t06SJJtdoBAABOJuDgMnLkSB04cEDZ2dkqKSnRgAEDtHbtWt8Ju3v37lVICBfkBQAATS/g67gEQ3N8DxxoqApPlS7IXieJ66oAQCCCfh0XAACAYCK4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANToEu4C2yhijY8e9wS4DTaDCwzwCQGtBcGkGxhiNyC9U0Z6DwS4FAIA2hUNFzeDYcS+hpQ1KjO+q8FBnsMsAgHaNPS7N7JOZKYpw8WbXFoSHOuVwOIJdBgC0awSXZhbhcirCxdMMAEBT4FARAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAajQou8+fPV0JCgsLCwpSUlKTNmzfX23fRokUaMmSIunbtqq5duyolJeWE/QEAAOoTcHBZsWKFMjMzlZOTo08//VT9+/dXamqq9u/fX2f/jRs36rbbbtOGDRtUWFiouLg4/eY3v9E333xzysW3NsYYVXiqVOHxBrsUAADaJIcxxgQyICkpSZdccomeffZZSVJ1dbXi4uJ0zz33KCsr66TjvV6vunbtqmeffVajR4+us09lZaUqKyt998vLyxUXF6fDhw8rMjIykHJbjDFGI/ILVbTnoF/71lmpinB1CFJVAAAET3l5uaKiopr0/TugPS4ej0dFRUVKSUn51wJCQpSSkqLCwsIGLaOiokLHjx/XaaedVm+f3NxcRUVF+W5xcXGBlBkUx457a4WWxPiuCg91BqkiAADanoB2BZSVlcnr9SomJsavPSYmRtu2bWvQMqZOnapevXr5hZ+fmzZtmjIzM333a/a42OKTmSmKcDkVHuqUw+EIdjkAALQZLXoM4/HHH9fy5cu1ceNGhYWF1dvP7XbL7Xa3YGVNK8Ll5PAQAADNIKB31+joaDmdTpWWlvq1l5aWqkePHicc+9RTT+nxxx/Xu+++q379+gVeKQAAaPcCOsfF5XJp4MCBKigo8LVVV1eroKBAycnJ9Y574oknNHv2bK1du1aJiYmNrxYAALRrAR/PyMzMVEZGhhITEzVo0CDl5eXp6NGjGjt2rCRp9OjRio2NVW5uriTp3//935Wdna1XXnlFCQkJKikpkSR16tRJnTp1asJNAQAAbV3AwWXkyJE6cOCAsrOzVVJSogEDBmjt2rW+E3b37t2rkJB/7ch57rnn5PF4NGLECL/l5OTk6OGHHz616gEAQLsS8HVcgqE5vgfe1Co8Vboge50krt0CAIDUCq7jAgAAEEwEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1ugQ7ALaivBQp7bOSvX9GwAAND2CSxNxOByKcPF0AgDQnDhUBAAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABYg+ACAACsQXABAADWILgAAABrEFwAAIA1CC4AAMAaBBcAAGANggsAALAGwQUAAFiD4AIAAKxBcAEAANYguAAAAGsQXAAAgDUILgAAwBoEFwAAYA2CCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgjUYFl/nz5yshIUFhYWFKSkrS5s2bT9j/tdde03nnnaewsDBddNFFWrNmTaOKBQAA7VvAwWXFihXKzMxUTk6OPv30U/Xv31+pqanav39/nf03bdqk2267TXfccYe2bNmi4cOHa/jw4fr73/9+ysUDAID2xWGMMYEMSEpK0iWXXKJnn31WklRdXa24uDjdc889ysrKqtV/5MiROnr0qN5++21f269+9SsNGDBA+fn5DVpneXm5oqKidPjwYUVGRgZSLgAACJLmeP/uEEhnj8ejoqIiTZs2zdcWEhKilJQUFRYW1jmmsLBQmZmZfm2pqalatWpVveuprKxUZWWl7/7hw4cl/fgEAAAAO9S8bwe4j+SEAgouZWVl8nq9iomJ8WuPiYnRtm3b6hxTUlJSZ/+SkpJ615Obm6tHHnmkVntcXFwg5QIAgFbgn//8p6KioppkWQEFl5Yybdo0v700hw4dUnx8vPbu3dtkG47GKS8vV1xcnL7++msO2wUZc9F6MBetC/PRehw+fFhnnnmmTjvttCZbZkDBJTo6Wk6nU6WlpX7tpaWl6tGjR51jevToEVB/SXK73XK73bXao6KieBG2EpGRkcxFK8FctB7MRevCfLQeISFNd/WVgJbkcrk0cOBAFRQU+Nqqq6tVUFCg5OTkOsckJyf79Zek9evX19sfAACgPgEfKsrMzFRGRoYSExM1aNAg5eXl6ejRoxo7dqwkafTo0YqNjVVubq4k6b777tPQoUM1d+5cXX/99Vq+fLk++eQTLVy4sGm3BAAAtHkBB5eRI0fqwIEDys7OVklJiQYMGKC1a9f6TsDdu3ev3y6hwYMH65VXXtHMmTM1ffp0/eIXv9CqVavUt2/fBq/T7XYrJyenzsNHaFnMRevBXLQezEXrwny0Hs0xFwFfxwUAACBY+K0iAABgDYILAACwBsEFAABYg+ACAACsQXABAADWaDXBZf78+UpISFBYWJiSkpK0efPmE/Z/7bXXdN555yksLEwXXXSR1qxZ00KVtn2BzMWiRYs0ZMgQde3aVV27dlVKSspJ5w4NF+jfRY3ly5fL4XBo+PDhzVtgOxLoXBw6dEgTJ05Uz5495Xa7de655/L/VBMJdC7y8vLUp08fhYeHKy4uTlOmTNEPP/zQQtW2Xe+//77S0tLUq1cvORyOE/54co2NGzfql7/8pdxut8455xy98MILga/YtALLly83LpfLLFmyxHz++efmzjvvNF26dDGlpaV19v/www+N0+k0TzzxhNm6dauZOXOmCQ0NNZ999lkLV972BDoXo0aNMvPnzzdbtmwxX3zxhRkzZoyJiooy//d//9fClbc9gc5FjV27dpnY2FgzZMgQc+ONN7ZMsW1coHNRWVlpEhMTzXXXXWc++OADs2vXLrNx40ZTXFzcwpW3PYHOxcsvv2zcbrd5+eWXza5du8y6detMz549zZQpU1q48rZnzZo1ZsaMGeb11183kswbb7xxwv47d+40ERERJjMz02zdutU888wzxul0mrVr1wa03lYRXAYNGmQmTpzou+/1ek2vXr1Mbm5unf1vvfVWc/311/u1JSUlmX/7t39r1jrbg0Dn4ueqqqpM586dzbJly5qrxHajMXNRVVVlBg8ebBYvXmwyMjIILk0k0Ll47rnnzNlnn208Hk9LldhuBDoXEydONL/+9a/92jIzM82ll17arHW2Nw0JLg8++KC58MIL/dpGjhxpUlNTA1pX0A8VeTweFRUVKSUlxdcWEhKilJQUFRYW1jmmsLDQr78kpaam1tsfDdOYufi5iooKHT9+vEl/CbQ9auxczJo1S927d9cdd9zREmW2C42Zi7feekvJycmaOHGiYmJi1LdvXz322GPyer0tVXab1Ji5GDx4sIqKinyHk3bu3Kk1a9bouuuua5Ga8S9N9d4d8CX/m1pZWZm8Xq/vJwNqxMTEaNu2bXWOKSkpqbN/SUlJs9XZHjRmLn5u6tSp6tWrV60XJwLTmLn44IMP9Pzzz6u4uLgFKmw/GjMXO3fu1Hvvvafbb79da9as0Y4dO3T33Xfr+PHjysnJaYmy26TGzMWoUaNUVlamyy67TMYYVVVVacKECZo+fXpLlIyfqO+9u7y8XMeOHVN4eHiDlhP0PS5oOx5//HEtX75cb7zxhsLCwoJdTrty5MgRpaena9GiRYqOjg52Oe1edXW1unfvroULF2rgwIEaOXKkZsyYofz8/GCX1u5s3LhRjz32mBYsWKBPP/1Ur7/+ulavXq3Zs2cHuzQ0UtD3uERHR8vpdKq0tNSvvbS0VD169KhzTI8ePQLqj4ZpzFzUeOqpp/T444/r3XffVb9+/ZqzzHYh0Ln46quvtHv3bqWlpfnaqqurJUkdOnTQ9u3b1bt37+Ytuo1qzN9Fz549FRoaKqfT6Ws7//zzVVJSIo/HI5fL1aw1t1WNmYuHHnpI6enpGjdunCTpoosu0tGjRzV+/HjNmDHD70eB0bzqe++OjIxs8N4WqRXscXG5XBo4cKAKCgp8bdXV1SooKFBycnKdY5KTk/36S9L69evr7Y+GacxcSNITTzyh2bNna+3atUpMTGyJUtu8QOfivPPO02effabi4mLfbdiwYbryyitVXFysuLi4liy/TWnM38Wll16qHTt2+MKjJH355Zfq2bMnoeUUNGYuKioqaoWTmkBp+I3hFtVk792BnTfcPJYvX27cbrd54YUXzNatW8348eNNly5dTElJiTHGmPT0dJOVleXr/+GHH5oOHTqYp556ynzxxRcmJyeHr0M3kUDn4vHHHzcul8usXLnS7Nu3z3c7cuRIsDahzQh0Ln6ObxU1nUDnYu/evaZz585m0qRJZvv27ebtt9823bt3N48++miwNqHNCHQucnJyTOfOnc2f/vQns3PnTvPOO++Y3r17m1tvvTVYm9BmHDlyxGzZssVs2bLFSDLz5s0zW7ZsMXv27DHGGJOVlWXS09N9/Wu+Dv3HP/7RfPHFF2b+/Pn2fh3aGGOeeeYZc+aZZxqXy2UGDRpkPvroI99jQ4cONRkZGX79X331VXPuuecal8tlLrzwQrN69eoWrrjtCmQu4uPjjaRat5ycnJYvvA0K9O/ipwguTSvQudi0aZNJSkoybrfbnH322WbOnDmmqqqqhatumwKZi+PHj5uHH37Y9O7d24SFhZm4uDhz9913m4MHD7Z84W3Mhg0b6vz/v+b5z8jIMEOHDq01ZsCAAcblcpmzzz7bLF26NOD1OoxhXxkAALBD0M9xAQAAaCiCCwAAsAbBBQAAWIPgAgAArEFwAQAA1iC4AAAAaxBcAACANQguAADAGgQXAABgDYILAACwBsEFAABY4/8Dn3GdXDevh2QAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.distributions.empirical_distribution import ECDF\n",
    "\n",
    "x = np.random.uniform(size=10)\n",
    "ecdf = ECDF(x)\n",
    "\n",
    "# get unique sorted values of x and corresponding values of y\n",
    "ux, uy = np.unique(ecdf.x, return_index=False, return_inverse=False, return_counts=False), ecdf.y[np.unique(ecdf.y, return_index=True, return_inverse=False, return_counts=False)[1]]\n",
    "\n",
    "# add last point (1,1) to the step function\n",
    "ux = np.concatenate((ux, [1]))\n",
    "uy = np.concatenate((uy, [1]))\n",
    "\n",
    "# create step function with horizontal lines\n",
    "plt.plot(ux[:-1], uy[:-1], drawstyle='steps-pre')\n",
    "plt.plot([ux[-2], ux[-1]], [uy[-2], uy[-2]], color='C0', linestyle='--')\n",
    "plt.plot([ux[-1], ux[-1]], [uy[-2], uy[-1]], color='C0', linestyle='--')\n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(0, 1)\n",
    "plt.title('ECDF for uniform distribution')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Nonparametric Bootstrap\n",
    "\n",
    "* Implementation of bootstrap is a simulation exercise. \n",
    "* In an i.i.d. environment, $n$ observations are drawn with equal weight and **with replacement** from the realized sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Variants of Bootstrap Schemes\n",
    "\n",
    "* Block bootstrap: preserve dependence structure\n",
    "  * dependent dataset such as time series\n",
    "  * clustering data or networks\n",
    "\n",
    "* parametric bootstrap\n",
    "  * In regressions we fix the regressors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* wild bootstrap: [[Davidson, 2010]](https://www.tandfonline.com/doi/abs/10.1198/jbes.2009.07221?journalCode=ubes20)\n",
    "  * Heteroskedascity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   no         x\n",
      "0   1 -0.243709\n",
      "1   2 -1.518069\n",
      "2   3  1.884239\n",
      "3   4  0.784840\n",
      "4   5 -1.621343\n",
      "5   6 -0.333652\n",
      "6   7 -0.594756\n",
      "7   8  0.066655\n",
      "8   9 -1.781558\n",
      "    no         x\n",
      "0    3  1.884239\n",
      "1    9 -1.781558\n",
      "2    2 -1.518069\n",
      "3    4  0.784840\n",
      "4    1 -0.243709\n",
      "5    4  0.784840\n",
      "6    5 -1.621343\n",
      "7    3  1.884239\n",
      "8    3  1.884239\n",
      "9    9 -1.781558\n",
      "10   4  0.784840\n",
      "11   5 -1.621343\n",
      "12   7 -0.594756\n",
      "13   9 -1.781558\n",
      "14   7 -0.594756\n",
      "15   3  1.884239\n",
      "16   9 -1.781558\n",
      "17   1 -0.243709\n",
      "18   2 -1.518069\n",
      "19   8  0.066655\n",
      "20   3  1.884239\n",
      "21   4  0.784840\n",
      "22   8  0.066655\n",
      "23   7 -0.594756\n",
      "24   6 -0.333652\n",
      "25   1 -0.243709\n",
      "26   4  0.784840\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "n = 9  # real sample size\n",
    "real_sample = np.random.normal(size=n)  # the real sample\n",
    "d0 = pd.DataFrame({'no': np.arange(1, n+1), 'x': real_sample})\n",
    "print(d0)\n",
    "\n",
    "# bootstrap\n",
    "boot_Rep = 3  # bootstrap 3 times\n",
    "d_boot = []  # save the bootstrap sample\n",
    "for b in range(boot_Rep):\n",
    "    boot_index = np.random.choice(n, n, replace=True)\n",
    "    d_boot.append(pd.DataFrame({'no': boot_index+1, 'x': real_sample[boot_index]}))\n",
    "\n",
    "d_boot = pd.concat(d_boot, ignore_index=True)  # combine the bootstrap samples\n",
    "print(d_boot)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bootstrap Estimation\n",
    "\n",
    "\n",
    "* R package [boot](http://cran.r-project.org/web/packages/boot/index.html) provides a general function `boot()`.\n",
    "* `ply`-family functions for repeated simulations. \n",
    "\n",
    "* Bootstrap is convenient. \n",
    "  * Analytic formula of the variance of an econometric estimator can be complex to derive or code up."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Example\n",
    "\n",
    "* One of the most popular estimators for a sample selection model is @heckman1977sample's two-step method\n",
    "\n",
    "* Outcome equation\n",
    "\n",
    "$$\n",
    "y_i = x_i \\beta + u_i\n",
    "$$\n",
    "\n",
    "and the selection equation be\n",
    "\n",
    "$$\n",
    "D_i = z_i \\gamma + v_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* To obtain a point estimator, we simply run a Probit in the selection model, predict the probability of participation, and then run an OLS of $y_i$ on $x_i$ and $\\lambda (\\hat{D}_i)$ in the outcome model, where $\\lambda(\\cdot)$ is the inverse Mill's ratio. \n",
    "\n",
    "* From Heckman (1979)'s original paper, the asymptotic variance expression of the two-step estimator is very complicated.\n",
    "\n",
    "* Instead of following the analytic formula, we can simply bootstrap the variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Below is the function for a single bootstrap. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Implementation is just a repeated evaluation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bootstrap Test\n",
    "\n",
    "* Bootstrap is particularly helpful in inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "* Test a hypothesis about the population mean. \n",
    "\n",
    "* Use $t$-statistic\n",
    "* Distribution of the sample is either\n",
    "  * normal\n",
    "  * zero-centered chi-square \n",
    "\n",
    "* We will show that the bootstrap test size is\n",
    "more precise than that of the asymptotic approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import t\n",
    "\n",
    "def T_stat(Y, mu):\n",
    "    n = len(Y)\n",
    "    return np.sqrt(n) * (np.mean(Y) - mu) / np.std(Y, ddof=1)\n",
    "\n",
    "def boot_test(Y, alpha, boot_Rep):\n",
    "    n = len(Y)\n",
    "    boot_T = np.zeros(boot_Rep)\n",
    "\n",
    "    for r in range(boot_Rep):\n",
    "        indices = np.random.choice(n, n, replace=True)\n",
    "        resampled_Y = Y[indices]\n",
    "        boot_T[r] = abs(T_stat(resampled_Y, np.mean(Y)))\n",
    "\n",
    "    boot_critical_value = np.quantile(boot_T, 1 - alpha)\n",
    "    return boot_critical_value\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bootstrap World\n",
    "\n",
    "* The null hypothesis must be imposed no matter the hypothesized parameter is true value or not.\n",
    "\n",
    "* Bootstrap $t$-statistic is\n",
    "\n",
    "$$\n",
    "T^{*}_{n} = \\frac{\\sqrt{n} (\\bar{X^{*}} - \\bar{X}) } { s^{*}  }.\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Design the Statistic\n",
    "\n",
    "\n",
    "* In the bootstrap world the **true** distribution is $F_n$\n",
    "* The bootstrap $t$-statistic is centered around $\\bar{X}$, the sample mean of $F_n$\n",
    "* If we wrongly center it around the population mean $\\theta$, then the test will have no power when the null hypothesis is false."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The following chuck of code report the rejection probability from three decision rules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0615 0.0995 0.0435]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import t, norm, chi2\n",
    "\n",
    "def T_stat(Y, mu):\n",
    "    n = len(Y)\n",
    "    return np.sqrt(n) * (np.mean(Y) - mu) / np.std(Y)\n",
    "\n",
    "def boot_test(Y, boot_Rep):\n",
    "    n = len(Y)\n",
    "    boot_T = np.zeros(boot_Rep)\n",
    "\n",
    "    for r in range(boot_Rep):\n",
    "        indices = np.random.choice(n, n, replace=True)\n",
    "        resampled_Y = Y[indices]\n",
    "        boot_T[r] = abs(T_stat(resampled_Y, np.mean(Y)))\n",
    "\n",
    "    boot_critical_value = np.quantile(boot_T, 1 - alpha)\n",
    "    return boot_critical_value\n",
    "\n",
    "def compare():\n",
    "    if distribution == \"normal\":\n",
    "        X = np.random.normal(size=n)\n",
    "    elif distribution == \"chisq\":\n",
    "        X = chi2.rvs(n-1, size=n) - n + 1\n",
    "\n",
    "    t_value_X = T_stat(X, mu)\n",
    "\n",
    "    exact = abs(t_value_X) > t.ppf(1 - alpha/2, n-1)\n",
    "    asym = abs(t_value_X) > norm.ppf(1 - alpha/2)\n",
    "    boot_rule = abs(t_value_X) > boot_test(X, boot_Rep)\n",
    "\n",
    "    return np.array([exact, asym, boot_rule])\n",
    "\n",
    "n = 10\n",
    "distribution = \"normal\"\n",
    "boot_Rep = 199\n",
    "MC_rep = 2000\n",
    "alpha = 0.05\n",
    "mu = 0\n",
    "\n",
    "res = np.zeros((MC_rep, 3))\n",
    "\n",
    "for i in range(MC_rep):\n",
    "    res[i,:] = compare()\n",
    "\n",
    "print(np.mean(res, axis=0))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Results\n",
    "\n",
    "* The program reports the empirical size.\n",
    "* Nominal size of the test is 5%. \n",
    "* Bootstrap test is more accurate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "* When the underlying distribution is a $\\chi^2$, the exact distribution is difficult \n",
    "to derive analytically. \n",
    "* However, we can still compare the asymptotic size with the bootstrap size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0625 0.0945 0.0435]\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "distribution = \"chisq\"\n",
    "boot_Rep = 199\n",
    "MC_rep = 2000\n",
    "alpha = 0.05\n",
    "mu = 0\n",
    "\n",
    "res = np.zeros((MC_rep, 3))\n",
    "\n",
    "for i in range(MC_rep):\n",
    "    res[i,:] = compare()\n",
    "\n",
    "print(np.mean(res, axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Here the \"exact test\" is no longer exact. \n",
    "* The asymptotic test works fairly reasonable\n",
    "* Bootstrap is closer to the nominal size 5%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Reading\n",
    "\n",
    "* [ISLR](https://www.statlearning.com/): Ch 5.2\n",
    "* [CASI](https://hastie.su.domains/CASI/): Ch 10.2-4"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "rise": {
   "enable_chalkboard": true,
   "scroll": true,
   "theme": "serif"
  },
  "vscode": {
   "interpreter": {
    "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}