{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "import estimagic as em"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap Monte Carlo Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this juypter notebook, we perform a Monte Carlo exercise to illustrate the importance of using the cluster robust variant of the bootstrap when data within clusters is correlated. \n",
    "\n",
    "The main idea is to repeatedly draw clustered samples, get both uniform and clustered bootstrap estimates in these samples, and then compare how often the true null hypothesis is rejected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Generating Process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The true data generating process is given by\n",
    "\n",
    "$$ logit(y_{i,g}) = \\beta_0 + \\beta_1 (x_{i,g}) + \\epsilon_{i,g}, $$\n",
    "\n",
    "where the independent variable $x_{i,g} = x_i + x_g$ and the noise term $\\epsilon_{i,g} = \\epsilon_i + \\epsilon_g$ each consist of an individual and a cluster term.\n",
    "\n",
    "In the simulations we perform below, we have $\\beta_0 = \\beta_1 =0$. $x_i$ and $x_g$ are drawn from a standard normal distribution, and $\\epsilon_i$ and $\\epsilon_g$ are drawn from a normal distribution with $\\mu_0$ and $\\sigma=0.5$. The value of $\\sigma$ is chosen to not blow up rejection rates in the independent case too much."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_clustered_data(nclusters, nobs_per_cluster, true_beta=0):\n",
    "    \"\"\"Create a bivariate clustered dataset with specified number of\n",
    "    clusters and number of observations per cluster that has a population\n",
    "    value of true_beta for the logit coefficient on the independent variable.\n",
    "\n",
    "    Args:\n",
    "        nclusters (int): Number of clusters.\n",
    "        nobs_per_cluster (int): Number of observations per cluster.\n",
    "        true_beta (int): The true logit coefficient on x.\n",
    "\n",
    "    Returns:\n",
    "        pd.DataFrame: Clustered dataset.\n",
    "    \"\"\"\n",
    "\n",
    "    x_cluster = np.random.normal(size=nclusters)\n",
    "    x_ind = np.random.normal(size=nobs_per_cluster * nclusters)\n",
    "    eps_cluster = np.random.normal(size=nclusters, scale=0.5)\n",
    "    eps_ind = np.random.normal(size=nobs_per_cluster * nclusters, scale=0.5)\n",
    "\n",
    "    y = []\n",
    "    x = []\n",
    "    cluster = []\n",
    "\n",
    "    for g in range(nclusters):\n",
    "\n",
    "        for i in range(nobs_per_cluster):\n",
    "\n",
    "            key = (i + 1) * (g + 1) - 1\n",
    "\n",
    "            arg = (\n",
    "                true_beta * (x_cluster[g] + x_ind[key]) + eps_ind[key] + eps_cluster[g]\n",
    "            )\n",
    "\n",
    "            y_prob = 1 / (1 + np.exp(-arg))\n",
    "            y.append(np.random.binomial(n=1, p=y_prob))\n",
    "            x.append(x_cluster[g] + x_ind[(i + 1) * (g + 1) - 1])\n",
    "            cluster.append(g)\n",
    "\n",
    "    y = np.array(y)\n",
    "    x = np.array(x)\n",
    "    cluster = np.array(cluster)\n",
    "\n",
    "    return pd.DataFrame({\"y\": y, \"x\": x, \"cluster\": cluster})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo Simulation Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes bootstrap t-values. As suggested my Cameron and Miller (2015), critical values are the 0.975 quantiles from a t distribution with `n_clusters` -1 degrees of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_t_values(data, sample_size=200, hyp_beta=0, cluster=False):\n",
    "    \"\"\"Get bootstrap t-values for testing the hypothesis that beta == hyp_beta.\n",
    "\n",
    "    Args:\n",
    "        data (pd.DataFrame): Original dataset.\n",
    "        sample_size (int): Number of bootstrap samples to draw.\n",
    "        hyp_beta (float): Hypothesised value of beta.\n",
    "        cluster (bool): Whether or not to cluster on the cluster column.\n",
    "\n",
    "    Returns:\n",
    "        float: T-Value of hypothesis.\n",
    "    \"\"\"\n",
    "\n",
    "    def logit_wrap(df):\n",
    "\n",
    "        y = df[\"y\"]\n",
    "        x = df[\"x\"]\n",
    "\n",
    "        result = sm.Logit(y, sm.add_constant(x)).fit(disp=0).params\n",
    "\n",
    "        return pd.Series(result, index=[\"constant\", \"x\"])\n",
    "\n",
    "    if cluster is False:\n",
    "\n",
    "        result = em.bootstrap(data=data, outcome=logit_wrap, n_draws=sample_size)\n",
    "        estimates = pd.DataFrame(result.outcomes())[\"x\"]\n",
    "\n",
    "    else:\n",
    "\n",
    "        result = em.bootstrap(\n",
    "            data=data,\n",
    "            outcome=logit_wrap,\n",
    "            n_draws=sample_size,\n",
    "            cluster_by=\"cluster\",\n",
    "        )\n",
    "        estimates = pd.DataFrame(result.outcomes())[\"x\"]\n",
    "\n",
    "    return (estimates.mean() - hyp_beta) / estimates.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def monte_carlo(nsim, nclusters, nobs_per_cluster, true_beta=0, n_cores=-1):\n",
    "    \"\"\"Run a Monte Carlo simulation for rejection rates and a logit data generating process.\n",
    "\n",
    "    Rejection rates are based on a t distribution with nclusters-1 degrees of freedom.\n",
    "\n",
    "    Args:\n",
    "        nsim (int): Number of Monte Carlo draws.\n",
    "        nclusters (int): Number of clusters in each generated dataset.\n",
    "        nobs_per_cluster (int) Number of observations per cluster.\n",
    "        true_beta (int): Population value of logit coefficient on x.\n",
    "        n_cores (int): Number of jobs for Parallelization.\n",
    "\n",
    "    Returns:\n",
    "        pd.DataFrame: DataFrame of average rejection rates.\n",
    "    \"\"\"\n",
    "\n",
    "    reject_independent = np.zeros(nsim)\n",
    "\n",
    "    reject_cluster = np.zeros(nsim)\n",
    "\n",
    "    def loop(i):\n",
    "\n",
    "        df = create_clustered_data(nclusters, nobs_per_cluster, true_beta)\n",
    "\n",
    "        return [get_t_values(df), get_t_values(df, cluster=True)]\n",
    "\n",
    "    t_value_array = np.array(\n",
    "        Parallel(n_jobs=n_cores)(delayed(loop)(i) for i in range(nsim))\n",
    "    )\n",
    "\n",
    "    crit = scipy.stats.t.ppf(0.975, nclusters - 1)\n",
    "\n",
    "    result = pd.DataFrame(np.abs(t_value_array) > crit, columns=[\"uniform\", \"cluster\"])\n",
    "\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we perform Monte Carlo simulations with the above functions. In each simulation, the sample size is 200, but the number of clusters varies across simulations. Be warned that the code below takes a long time to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(505)\n",
    "\n",
    "results_list = []\n",
    "\n",
    "for g, k in [[20, 50], [50, 20], [100, 10], [200, 5], [500, 2]]:\n",
    "\n",
    "    results_list.append(monte_carlo(nsim=500, nclusters=g, nobs_per_cluster=k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           uniform  cluster\n",
      "nclusters                  \n",
      "20           0.146    0.070\n",
      "50           0.080    0.050\n",
      "100          0.088    0.078\n",
      "200          0.066    0.054\n",
      "500          0.046    0.046\n"
     ]
    }
   ],
   "source": [
    "mean_rejection_data = pd.DataFrame([x.mean() for x in results_list])\n",
    "mean_rejection_data[\"nclusters\"] = [20, 50, 100, 200, 500]\n",
    "mean_rejection_data.set_index(\"nclusters\", inplace=True)\n",
    "\n",
    "print(mean_rejection_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'Comparison of Rejection Rates')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEmCAYAAACTYry7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABFTUlEQVR4nO3dd3hUZfbA8e8hdBBEjKK0ACYqNUixoYBYENvqih1ldWXZtbcVxd4LinVV1q641l3Wn70siF1AUQEFQUAQpIpIb+f3x7lDboZJZpJMSTLn8zzzZObW994k98zbRVVxzjnnSlMj0wlwzjlX+XmwcM45F5cHC+ecc3F5sHDOOReXBwvnnHNxebBwzjkXlweLSkxEjhOR/4nIChFZLyIzROQmEdkx02lLFhHJExEVkSMznZayEJE9ReRDEVkdpD+vhO3mBOtVRDaIyA8icruINCjHOVVEzq1w4osfs0BErhOR7aOWDw7O1zCZ5yslHYND90lFZKmIjBWRA8txrJjX5CrGg0UlJSJ3AS8BPwKDgEOBkcBRwD8zmLRkWwjsC3yU6YSU0Z3A9sDRWPoXlrLtc8E2BwNPAxcB95bjnPtifxPJVABci11L2OvB+dYk+XzxHBSc989AA+BtEckv4zFKuiZXATUznQC3LRE5CrgYOEtVHw+t+kBERmGBo8oTkbqqug74LNNpKYc9gFdV9f0Etl2oqpFrHC8izYHBIjJEVbckesLQMVJOVZcAS9J1vpAJqroKQES+AOYDA4FbMpAWF+I5i8rpIuDLqEABgKpuVtU3I59FZEcReUpElonIGhEZJyLdw/sERSEjRGSYiCwUkd9E5C4xA0Rkqoj8LiJjRKRJaL8+QZHAoSLyWlDk8pOIDI06/r4i8qqILAi2mSwip0ZtEylm6BmkcS1wWaxiKBE5WkQmBcf6VUQ+F5HeofX1ReQ+EflFRNaJyAQROTTqfONE5GUROUVEZorIShF5U0RaxLv5IlIoIu8H9/NXERktIjsH6/JERIF2wEVB2sfFO2aUr4E6QG7onDuIyCMisii4pk9EZO+odG1TDCUix4jIxGCfX0TkDhGpFbVNZxH5P7HizFUi8oWIHCIifYD/CzabHRx/TrDPNsVQZfxbu0hE5gf37/nyFAmp6gIsYLUMHX+P4HjzgjRMFZELRaRGsL7EawrWtwr2Xx7s/7aI7B51DVcEfzPrgt/HWyLSrKzpr248Z1HJBP/o+wF3JbjLGGA34FJgKXAZMFZEuqrqzNB2JwFfAH8CugE3YV8WDgSuBuoBDwC3AsWCAfAY8AxwP3Ac8JCIzFfV14L1rYGPgYeBdcD+wBMiskVV/xV1rH8BDwHXAytiXH874GWsmOYyoG6Q3h1Cm/0TK/65EpgJnA28LiJ9VTVcnLU3sCtwSXB99wKjgAHR5w2dPxcYB3wHnAI0BG4D3g0ejJFis/8A/wvuycqSjleCVsDv2O8LEakDvIcVm1wGLAb+CrwnIvmq+ksJaT0Bu5+PYPeiHfb7q4H9PSAie2C/m+nY73UZ0B17AL8cbDcC+70uBNaXku4xJPa3dgLwDTAEaAHcjeUM/hb/1hS7vobY7312aHHz4FpGY/ewEPtbqhdc+5clXZOI7IAVdy4L7sUaYBh2nwtUda2InI7dy8uBqUBTrGiszHVM1Y6q+qsSvYBmgAJ/SWDb/sG2vUPLGmDfxh4JLZuDPVRzQsu+ADYBbULL7gAWhT73CY4/Kuq87wKflZAmwb6EPAL8L7R8cHCsC6K2zwuWHxl8Ph5YVso17wlsAc4ILasBTAHeDi0bB/wGNAktuzA4V71Sjn8bFsQahZb1DPY7OeqejkjgdzQHC/w1gfrB7+xX4PLQNmcBG4D80LKawCzgztAyBc4N3ee5wBNR5zsTWAs0DT7/CyvKiXnNwJHBcfOilkd+Xw3L8bc2C6gZWnYP8Euc+xQ5X+Pg2psDzwI/Ablx/tauBH5M4JpuxALFDqFlTYK/k3OCzw8Ar1Tkf7i6vrwYqvJKZITHnsASVf1g606qq4HXgF5R245T1c2hzzOBOao6O2pZrojUjtr3P1Gf/w10E5EcABFpEhQLzQU2Bq8hWEVjtNfjXNO3QOOguONQ2bbVUA/sIbG1olet3P8ltr3mCar6a+jztOBn81LO3xN4R1W35hZU9QvsIRh9/ERdjN2T1cCbwFhVvT20/mBgElZsUlNEIjn+D7BcQCwFWA7lxcg+wX7/w3JjHYPtDgJeUNW15Ux7RFn+1saq6qbQ52nATjH+rmJZgd2r+cAfgePV6k8Aq+cSketFZCaWY9gI3Ay0Cd23khyMfdFZGbpfv2P3PnKfJwMDgnP0jPyNO6+zqIyWYf8ErRLYdhdgUYzliyhebAPbFvlsKGGZANH/1ItjfK4JRJrwPgmciLUQOhR7oD+OPbRipa1EqjodOAZoC7wBLBWR54LiIbBrXqWq0a10FgH1gyKdiBVR22wIfsZKV0RZ7mminsXuSR/gCeBYEflraP2OwD4UBdrI60+EyuujRO79G1H7RIJ/ZL+mlN5SK1EV/VuL9XcVy4FYYDoNWA48H/WF4XasmClSnNgDK1KF0n+vYPfsRLa9z30pul+PYzmVE4DPgUUicqMHDa+zqHRUdaOIfAwcBlwVZ/OFwE4xlu+M/aMlS/Q5dsKKsJaKSF3gCKx45OHIBpEKxxji5phU9XWsDqJxcOx7sLqBk7Brbigi9aMCxs7AGlUtrcw9EaXd00nlPOYiVZ0YvP9ARFoDN4jI08G38+XARKyeIlpJ1xP5/Q4BvoqxPhI0lmEP+opK19/aV2qtoSaIyE/AeOBcLEiAtYy6X1XviOwgIkckeOzlwKtYcVS032FrLnUkMFJEWgKnYjmXn7E6uazlOYvK6R6gu4icEb1CRGqISP/g4+dY9v7A0Pr62AM2mf0Wjo3xeVJQrFUHyCH0UBOR7bAK6ApR1d9U9TmsGKx9sHgCFnCOD51Pgs/JuObPgcOCa4gcvwdWt5Kse3oF9i33rODz+1jF8U+qOjHq9W0Jx5iOPcDyYuwzUVWXhY59QhDUY0kktwXp+1vbSlU/xHJOF4ZyjPUo/reWg32JCCvpmt4HOgBTY9yv6THOP09Vb8OKZ9tHr882nrOohFT1/0TkbuAxEdkf+C+wCmvbPxQrP39LVd8OciEviMgw7Fvkpdg/1J1JTNLhInIzVoZ+HHAIVlSEqv4mIhOAa0RkJVb5PAyrNGxU1hOJyF+w1kZvAQuAfOzb5NPB+b4TkX8BD4hII4paQ+1B7G/mZXV3cJy3ReR2ilpDfQu8koTjo6pfiMi7WNPbB7FrGwqME5ERWEfMplhxzC+qOjLGMbaIyCXAM8F9eBN7SLYF/oCV9a/BWgpNwPp33IX9jXTFGhE8jgUdgL+IyPNY7mybAJXGv7VoN2OtuU7DWuW9C5wT1FksB87BvrCElXRNdwfH+Z+I3I8F252B3sBHqvovEXkkOO5n2N9wX+xv8PLUXWIVkekadn+V/MIq+MZif7QbgBlYk8BmoW1ysYfNr1grmA+AHlHHmUNUyx2snmFi1LLBFG8B0yf4fBj2MFqDVTz+LWq/3bCK1dVY65W/A9cBS0s6dmh5HsVbQ+2LVYIvwJrhzsaKIOqE9qmPFUstwr5lTgQOizruOODlqGWR6+kY5753Da5nDVb+/hywc7x7WsKxYm6Hlc1vbWGFtQK6F5gX/K7nYw0J9g/towStdkLLDgc+DO79SqyC9iaKt0bqjH1D/z14fQ70C62/BGtZtQlr9BDz90X5/9Zi/u4T3Sb4XXyH1XvsjOU0Vwa//zuwLwvRad3mmoLlu2L1RpG/nTlYnVKHUDo+xgLGGqwJ8FmZfhZUhpcEN8i5bQQdnMYCnVR1SmZTk92CYrGVwCm6bd8V51LOi6Gcq+SCjnWnY9+ev8xwclyW8mDhXOV3K1Z/cZnGqIh1Lh28GMo551xc3nTWOedcXB4snHPOxeXBwjnnXFweLJxzzsXlwcI551xcHiycc87F5cHCOedcXB4snHPOxeXBwjnnXFweLJxzzsXlwcI551xcHiycc87FVW1Hnd1xxx01Ly8v08lwzrkqY9KkSUtVNTfWumobLPLy8pg4cWKmk+Gcc1WGiMwtaZ0XQznnnIvLg4Vzzrm4PFg455yLq9rWWTjn4tu4cSPz589n3bp1mU6KS6O6devSokULatWqlfA+Hiycy2Lz589nu+22Iy8vDxHJdHJcGqgqy5YtY/78+bRp0ybh/bwYKmT0aMjLgxo17Ofo0ZlOkXOptW7dOpo2beqBIouICE2bNi1zbtJzFoHRo2HIEFizxj7PnWufAU49NXPpci7VPFBkn/L8zj1nERg+vChQRKxZY8udcy7bebAI/PRT2ZY75ypuzpw5dOzYsdiy6667jhEjRpS638SJEzn//PMBWL9+PQcffDCFhYW88MILKUsrQE5ODoWFhXTp0oW99tqLTz75pFzHmTNnDs8991yp26xYsYJ//OMf5Tp+KniwCLRqVbblzmWjylKv1717d+677z4AvvrqKzZu3MjkyZM58cQTE9p/8+bN5TpvvXr1mDx5Ml9//TW33norV1xxRbmOU9FgUd70V4QHi8DNN0P9+sWX1a9vy51zRfV6c+eCalG9XioDRp8+fbj88svp2bMnBQUFfPjhhwCMGzeOI488ksWLF3PaaacxefJkCgsLmTVrFu+//z5du3alU6dOnHnmmaxfvx6wIYBuuOEGevXqxUsvvUReXh5XXnkl++67L927d+fLL7/ksMMOo127djz88MNx07Zy5UqaNGkCWAujyy67jI4dO9KpU6etOZySlg8bNowPP/yQwsJCRo4cydSpU+nZsyeFhYV07tyZH374gWHDhjFr1iwKCwu57LLLGDduHH379uWUU06hU6dOAPzhD3+gW7dudOjQgVGjRm1NW8OGDbnkkkvYa6+96NevH0uWLKn4L0NVq+WrW7duWlbPPqvauLEqqLZsaZ+dq86mTZu29f0FF6j27l3yq04d+9+IftWpU/I+F1xQ+vlnz56tHTp0KLbs2muv1TvvvFNVVXv37q0XX3yxqqq+/vrr2q9fP1VVHTt2rB5xxBHbvF+7dq22aNFCp0+frqqqgwYN0pEjR6qqauvWrfX222/fep7WrVvrP/7xD1VVvfDCC7VTp066cuVKXbx4sebm5sZMb40aNbRLly66++67a6NGjXTixImqqvryyy/rwQcfrJs2bdJffvlFW7ZsqQsWLChxeTjNqqrnnnuuPhs8cNavX69r1qzZ5t6MHTtW69evrz/++OPWZcuWLVNV1TVr1miHDh106dKlqqoKbD3e9ddfr+ecc8421xL+3UcAE7WEZ6rnLEJOPRUeesjev/mmt4JyLiz4gp7w8kSU1ConvPy4444DoFu3bsyZM6fU402fPp02bdpQUFAAwBlnnMH48eO3ro8upjr66KMB6NSpE3vvvTfbbbcdubm51K1blxUrVmxz/Egx1Pfff89bb73F6aefjqry0UcfcfLJJ5OTk8POO+9M7969mTBhQonLo+27777ccsst3H777cydO5d69erFvL6ePXsW6xtx33330aVLF/bZZx/mzZvHDz/8AECNGjW2Xutpp53GRx99VOp9S4Q3nY2Sn28/Z8yADh0ymxbn0umee0pfn5dnRU/RWreGcePKd86mTZvy66+/Flu2fPnyYg/EOnXqAFa5vGnTplKPZ1+OS9agQYNinyPHrlGjxtb3kc/xzrXvvvuydOlSlixZUuJ546Un4pRTTmHvvffm9ddf57DDDuPRRx+lbdu2paZ/3LhxvPfee3z66afUr1+fPn36lNh3IhnNoz1nESUSLIIA7ZwLpKJer2HDhuyyyy68//77gAWKt956i169epXreHvssQdz5sxh5syZADzzzDP07t27/Aksxffff8/mzZtp2rQpBx54IC+88AKbN29myZIljB8/np49e5a4fLvttuP333/feqwff/yRtm3bcv7553P00UfzzTffbLNNtN9++40mTZpQv359vv/+ez777LOt67Zs2cLLL78MwHPPPVfu+xnmOYsojRvDTjt5sHAuWqRYdvhwa1LeqpUFiooW1z799NOcc845XHLJJQBce+21tGvXrlzHqlu3Lk888QQDBw5k06ZN9OjRg6FDh1YsgSFr166lsLAQsFzDU089RU5ODsceeyyffvopXbp0QUS44447aNasWYnLmzZtSs2aNenSpQuDBw9m3bp1PPvss9SqVYtmzZpxzTXXsMMOO7D//vvTsWNHDj/8cI444ohiaenfvz8PP/wwnTt3Zvfdd2efffbZuq5BgwZMnTqVbt260bhx46Q0KZZEs0lVTffu3bW8kx/16gU5OfDBB0lOlHOVzHfffceee+6Z6WS4JGvYsCGrVq0qdZtYv3sRmaSq3WNtn7ZiKBHpLyLTRWSmiAyLsX4PEflURNaLyKUx1ueIyFci8lqq05qf7zkL55wLS0uwEJEc4EHgcKA9cLKItI/abDlwPlBS180LgO9SlsiQggJYuBDiBGbnnKuU4uUqyiNdOYuewExV/VFVNwDPA8eEN1DVxao6AdgYvbOItACOAB5NR2K9kts554pLV7BoDswLfZ4fLEvUPcDfgS2lbSQiQ0RkoohMrEiPRQ8WzjlXXLqCRaxGvgnVrIvIkcBiVZ0Ub1tVHaWq3VW1e25ublnTuNVuu9lPDxbOOWfSFSzmAy1Dn1sACxLcd3/gaBGZgxVfHSQizyY3ecU1aADNm1vHPOecc+kLFhOAfBFpIyK1gZOAVxPZUVWvUNUWqpoX7Pc/VT0tdUk13iLKufT45ZdfOOmkk2jXrh3t27dnwIABzJgxI+bw5Yl68sknWbAg0e+jJR8jNzeXwsJCOnTowPHHH8+a6ElvEjRmzBimTZtW6jbjxo0r95Dn6ZCWYKGqm4BzgbexFk0vqupUERkqIkMBRKSZiMwHLgauEpH5ItIoHemLpaDAg4Vz20jyGOWqyrHHHkufPn2YNWsW06ZN45ZbbmHRokUVOm55gkWs4T1OPPFEJk+ezNSpU6ldu3a5O7dVNFjEG3okHdLWz0JV31DVAlVtp6o3B8seVtWHg/e/BDmIRqq6ffB+ZdQxxqnqkelIb34+LF0KUcPWOJe9UjBG+dixY6lVq1axXtaFhYUccMABxbZ78sknOffcc7d+PvLIIxk3bhybN29m8ODBW4cAHzlyJC+//DITJ07k1FNPpbCwkLVr1zJp0iR69+5Nt27dOOyww1i4cCFgQ6BfeeWV9O7dm3vvvbfEdG7atInVq1dvHZJ87ty59OvXj86dO9OvXz9+CmZJi7X8k08+4dVXX+Wyyy7bOoz6fffdR/v27encuTMnnXQSc+bM4eGHH2bkyJEUFhby4YcfMnjwYC6++GL69u3L5ZdfzhdffMF+++1H165d2W+//Zg+ffrWe3PMMcfQv39/dt99d66//vpy/z5K48N9lCDcIqpnz8ymxbm0uPBCmDy55PWffbbtELNr1sBZZ8E//xl7n8LCUkconDJlCt26dStjQotMnjyZn3/+mSlTpgA2YdD222/PAw88wIgRI+jevTsbN27kvPPO47///S+5ubm88MILDB8+nMcff3zrPh+UMFzDCy+8wEcffcTChQspKCjgqKOOAuDcc8/l9NNP54wzzuDxxx/n/PPPZ8yYMSUuP/rooznyyCM5/vjjAbjtttuYPXs2derU2ZrmoUOH0rBhQy691PokP/bYY8yYMYP33nuPnJwcVq5cyfjx46lZsybvvfceV155Ja+88goAX3zxBVOmTKF+/fr06NGDI444gu7dY3bELjcfSLAEwQjHXhTlXEQqxiivoLZt2/Ljjz9y3nnn8dZbb9Go0bYl19OnT2fKlCkccsghFBYWctNNNzF//vyt60ubXS9SDPXLL7/QqVMn7rzzTgA+/fRTTjnlFAAGDRq0dQjwkpZH69y5M6eeeirPPvssNWuW/J194MCB5OTkADZw4MCBA+nYsSMXXXQRU6dO3brdIYccQtOmTalXrx7HHXdcUoYkj+Y5ixK0bWvFst4iymWNDIxR3qFDh62jo5amZs2abNlS1M0qMhR3kyZN+Prrr3n77bd58MEHefHFF7fmGCJUlQ4dOvDpp5/GPHb0sOWxiAhHHXUU999/P8OGbTNaUULzcoS9/vrrjB8/nldffZUbb7yx2IO/pLRdffXV9O3bl//85z/MmTOHPn36lHieZAxJHs1zFiWoU8dG1fSchXOBFIxRftBBB7F+/Xr+GSrGmjBhwjbFQnl5eUyePJktW7Ywb948vvjiCwCWLl3Kli1b+OMf/8iNN97Il19+CVBseO/dd9+dJUuWbA0WGzduLPHhXJqPPvpo62i4++23H88//zwAo0eP3joEeEnLw+mJXEPfvn254447WLFiBatWrUpoSPLmza0v85NPPlls3bvvvsvy5ctZu3YtY8aMYf/99y/z9cVV0hR6Vf1VnmlVox16qGr37hU+jHOVVqypNUv17LOqrVuritjPJMw9/PPPP+vAgQO1bdu22r59ex0wYIDOmDGj2LSiW7Zs0VNOOUXbt2+vJ5xwgvbu3VvHjh2rkydP1q5du2qXLl20S5cu+sYbb6iqTXNaUFCgXbp00TVr1uhXX32lBxxwgHbu3Fnbt2+vo0aNUlWbtnXChAkx0/XEE0/ojjvuqF26dNFOnTrp4YcfrosWLVJVmw62b9++2qlTJz3ooIN07ty5pS7/6KOPdM8999TCwkL9/vvvdf/999eOHTtqhw4d9NZbb1VV1enTp2unTp20S5cuOn78eD3jjDP0pZde2pqeTz75RPPz83W//fbTq666Slu3br01nQMHDtQBAwZoQUGBXnfddQnd97JOq+pDlJfi3HPhmWdgxQpIQa7OuYzzIcqrvieffJKJEyfywAMPlGm/SjtEeVWUnw8rV0IFhplyzrlqwYNFKbxFlHOushs8eHCZcxXl4cGiFJG+Ft4iylVn1bUo2pWsPL9zDxalyMuDmjU9Z+Gqr7p167Js2TIPGFlEVVm2bBl169Yt037ez6IUNWtafwsPFq66atGiBfPnz6ci87+4qqdu3bq0aNGiTPt4sIgjP9+LoVz1VatWLdq0aZPpZLgqwIuh4sjPh5kzbdw055zLVh4s4igosLHSKjg0vnPOVWkeLOLwFlHOOefBIq7wUOXOOZetPFjE0bKlDSrowcI5l808WMRRowbstpsXQznnspsHiwTk53vOwjmX3TxYJKCgAGbNgs2bM50S55zLDA8WCcjPhw0bIJiT3Tnnso4HiwR4iyjnXLbzYJEAH6rcOZftPFgkoFkzaNjQW0Q557KXB4sEiFjzWc9ZOOeylQeLBBUUeLBwzmUvDxYJys+H2bNh48ZMp8Q559LPg0WC8vOtn8Xs2ZlOiXPOpZ8HiwR5iyjnXDbzYJEgH6rcOZfN0hYsRKS/iEwXkZkiMizG+j1E5FMRWS8il4aWtxSRsSLynYhMFZEL0pXmsKZNYfvtPWfhnMtOaZmDW0RygAeBQ4D5wAQReVVVp4U2Ww6cD/whavdNwCWq+qWIbAdMEpF3o/ZNORFvEeWcy17pyln0BGaq6o+qugF4HjgmvIGqLlbVCcDGqOULVfXL4P3vwHdA8/Qku7j8fC+Gcs5lp3QFi+bAvNDn+ZTjgS8ieUBX4PPkJKts8vNh3jxYty4TZ3fOucxJV7CQGMu0TAcQaQi8AlyoqitL2GaIiEwUkYlLliwpRzJLV1AAqjZcuXPOZZN0BYv5QMvQ5xbAgkR3FpFaWKAYrar/Lmk7VR2lqt1VtXtubm65E1sSbxHlnMtW6QoWE4B8EWkjIrWBk4BXE9lRRAR4DPhOVe9OYRrj8qHKnXPZKi2toVR1k4icC7wN5ACPq+pUERkarH9YRJoBE4FGwBYRuRBoD3QGBgHfisjk4JBXquob6Uh7WOPGsNNOHiycc9knLcECIHi4vxG17OHQ+1+w4qloHxG7ziMjvEWUcy4beQ/uMsrP95yFcy77eLAoo4ICWLgQVq3KdEqccy59yhQsgqE39klVYqoCr+R2zmWjhIKFiLQSkY+B74H3gmXHi8ijqUxcZeTBwjmXjRLNWTwCvA5sR9FwHO9iYz1lld12s58eLJxz2STR1lA9gSNUdYuIKICq/iYijVOXtMqpQQNo3txbRDnnskuiOYtFwG7hBSLSHvgp6SmqArxFlHMu2yQaLEYAr4nIn4CaInIy8AJwe8pSVon5UOXOuWyTUDGUqj4uIsuBIdjosacDV6vqmBSmrdLKz4elS+HXX6FJk0ynxjnnUi+hYCEieweBYUzU8p6q+kUK0lWphVtE9eyZ2bQ451w6JFoM9W4Jy99KVkKqkoIC++lFUc65bFFqzkJEamDjMkkw+mt4jKZ22JSnWadtW6hRw1tEOeeyR7xiqE0UTVIUHRi2ADcnPUVVQJ060KqV5yycc9kjXrBog+UmPgAODC1XYImqrk1Vwio7bxHlnMsmpQYLVZ0bvG2dhrRUKfn58MwzNs2qVJoB1J1zLjUSns9CRI4GegM7Eqq7UNXTU5CuSi8/H1auhCVLbEIk55yrzhIdSPBabHyoGsBAYBlwGLAiZSmr5LxFlHMumyTadPZM4BBVvQjYEPw8CshLVcIqu0hfC28R5ZzLBokGi+1VdUrwfoOI1Ao64/VOUboqvbw8qFnTcxbOueyQaJ3FLBHpoKpTgSnAX0XkV+DX1CWtcqtZ0/pbeLBwzmWDRIPFVUDT4P0w4DmgIfC3VCSqqsjP92Io51x2iBssgl7c64DPAILip91K3SlL5OfD2LHefNY5V/3FrbNQ1S3Af1V1QxrSU6UUFMCaNbBgQaZT4pxzqZVoBfd4EdknpSmpgrxFlHMuWyRaZzEXeFNE/ovNZxEZLwpVvSYVCasKwkOV9+2b2bQ451wqJRos6lE0l0WL0HLddtPs0bKlDSroLaKcc9VdojPl/SnVCamKatSA3XbzYijnXPWXaJ2FK0F+vucsnHPVnweLCsrPh1mzYPPmTKfEOedSx4NFBRUUwIYN8NNPmU6Jc86ljgeLCgq3iHLOueoq4WAhIruLyAkicmb4VYb9+4vIdBGZKSLDYqzfQ0Q+FZH1InJpWfbNJA8WzrlskFBrKBG5ErgG+BpYE1qlwOMJ7J8DPAgcAswHJojIq6o6LbTZcuB84A/l2DdjdtkFGjTwFlHOueot0X4WFwI9VfWbcp6nJzBTVX8EEJHngWOArQ98VV0MLBaRI8q6byaJeIso51z1l2gx1Frg+wqcpznW8ztifrAsqfuKyBARmSgiE5csWVKuhJaHBwvnXHWXaLC4GrhfRHYRkRrhV4L7xxqTNdHe3wnvq6qjVLW7qnbPzc1N8PAVV1AAs2fDxo1pO6VzzqVVog/7J4GzsW/1G4PXpuBnIuYDLUOfWwCJjtVakX3TIj/f+lnMnp3plDjnXGokWmfRpoLnmQDki0gb4GfgJOCUNOybFuEWUQUFmU2Lc86lQqJjQ82FrRMh7QwsCua5SIiqbhKRc4G3gRzgcVWdKiJDg/UPi0gzYCLQCNgiIhcC7VV1Zax9E77CNIgEiBkz4Ijo6nnnnKsGEm062wh4APtWXxPYGLRKOl9Vf0vkGKr6BvBG1LKHQ+9/ofiItqXuW5k0bQrbb++V3M656ivROov7gAZAR2y48k5A/WB51vPms8656i7ROov+QFtVjXTImyEifwJmpSZZVU9BAXz4YaZT4ZxzqZFozmIdEN0WdUdgfXKTU3Xl58O8ebBuXaZT4pxzyZdosHgUeFdEhorI4UHF9NvAqNQlrWrJzwdVG67cOeeqm0SLoW7G+jacAuwavL+DBMaFyhbhFlEdOmQ2Lc45l2yJNp2NDBjowaEEPvqsc646KzFYiMggVX0meF/iUOSq6gEEaNwYcnM9WDjnqqfSchYnA88E7weVsE1CQ5Rni4ICH6rcOVc9lRgsVHVA6H3f9CSnasvPh7ffznQqnHMu+RJqDSUiX5WwfGJyk1O15efDwoWwalWmU+Kcc8mVaNPZ3aIXiIgAbZObnKot0iLK6y2cc9VNqa2hROTp4G3t0PuIPKBSDeiXaeEWUV27ZjYtzjmXTPGazs4q4b0CHwMvJT1FVdhuQf7LcxbOueqm1GChqtcDiMhnqupVt3E0aADNm3uLKOdc9ZNonUVXEekRXiAiPUXk7ylIU5Xmo88656qjRIPFBcC0qGXTgAuTmppqoLIEi9GjIS8PatSwn6NHZzpFzrmqLNGxoWqz7XzbG4C6yU1O1VdQAEuXwq+/QpMmmUnD6NEwZAisCQaUnzvXPgOcempm0uScq9oSzVlMAv4WtWwo8GVyk1P1VYYxoq64oihQRKxZA8OHZyY9zrmqL9GcxUXYEOWDsFZRu2FzcR+SqoRVVeFg0bNnes+9Zg088ojNqxHLTz+lNz3Oueoj0VFnp4pIAXAk0BL4N/Caqnpf5Sjt2tk0q+lsERUJErffDosWQZ06sD7GtFQ77pi+NDnnqpdEi6EIAsPHwIeq+rwHitjq1IHWrdNTDLVmDdx9N7RtCxdfDB07wvjx8NhjUL9+8W1FYMkSq7v4/ffUp805V70kOjZUKxH5GPgeeC9YdryIPJrKxFVVqW4RFQkSbdrAJZcUBYn33oMDDrBK7FGjLGiJ2M8nnoDLL4dHH4UuXWx755xLVKI5i0eA14HtKGoV9S5eZxFTZKhy1eQed/VquOuuoiDRuTN8+GFRkAg79VSYMwe2bLGfZ5wBt91mQUIE+vSBSy/1OcOdc4lJNFj0BG5T1S3YUB+o6m9A41QlrCrLz4eVK63YJxlWr4YRIyxIXHqp5Qw++gjefRd69SrbsXr1gq+/hr/8xQJPt27wpbdpc87FkWiwWETUyLMi0h7w9jUxJKv5bDhIXHYZFBZakHjnHdh///Ift2FDeOghePNNWLEC9t4bbrwRNm2qWHqdc9VXosFiBPCaiPwJqCkiJwMvALenLGVVWGSo8vK2iFq9Gu68syhIdO0KH39c8SARrX9/+PZbGDgQrrkG9tsPvv8+ecd3zlUfCQWLYJ7tvwMDgXnAGcDVquqDSMSQlwc1a5Y9Z7FqFdxxhwWJv/+9KEi8/bY9yFNhhx3guefghRdg1iw75333WV2Hc85FlKXp7BhVHaCqHVS1v6qOSWG6qrSaNe2Bn2iwCAeJyy+HvfaCTz5JbZCIdsIJMGUKHHQQXHABHHywd+JzzhUpsVOeiAxS1WeC92eWcowNwGxV/TjZiavKIi2iSrNqFTz4oNVLLF1qxULXXgv77JOeNEbbZRd47TXrp3HRRdCpk+UyTj/dWlA557JXaT24TwaeCd4PKmW7GkBbERmrqqcnLWVVXH4+jB1rzWejH7SRIHHnnbBsWeaDRJgI/PnP0K+fNbcdPBj+8x/rt7HTTplOnXMuU0oshlLVAaH3fUt59QZ2B44p7UQi0l9EpovITBEZFmO9iMh9wfpvRGSv0LqLRGSqiEwRkX+JSKUf7TY/3zrPLVhQtOz3362vQ14eDBtmY0d99pm1SqoMgSKsTRsLdiNGwFtvQYcOFjScc9kp4ToLEWkqIoNE5LLg864i0gJAVdcQ1bQ2at8c4EHgcKA9cHLQ9DbscCA/eA0BHgr2bQ6cD3RX1Y5ADnBSounOlMhgfi1bQqtWVieQl2cjwkaCxBtvWLPVyionxzr/TZpk13DccZbbWLEi0ylzzqVbosN99AamA6cC1wSL8wke6ACqWloXtJ7ATFX9UVU3AM+zbU7kGOBpNZ8B24vILsG6mkA9EakJ1AcWUImNHg333mvvVS1wvPSSBY7PP6/8QSJahw4W3K65xq6tUyfrNe6cyx6J5izuAU5U1f5ApOvW51gQSERzrMltxPxgWdxtVPVnrJ/HT8BC4DdVfSfB82bE8OGwdu22y1esSP+w5clSqxZcf7210mrQAA45BM47b9t5M5xz1VOiwSJPVd8P3kdGPNpA4vNhxGpLEz1yUsxtRKQJlutoA+wKNBCR02KeRGSIiEwUkYlLkjXWRjmU1OS0OjRF7dnThge54AJ44AHrl/H555lOlXMu1RINFtNE5LCoZQcD3ya4/3xsHoyIFmxblFTSNgdjTXOXqOpGbC6NmL0PVHWUqnZX1e65ubkJJi35WrUq2/KUSOEk3PXrwz33wPvv20CE++0HV10FGzYk7RTOuUom0WBxCTBaRJ7C6g4eAZ4ELktw/wlAvoi0EZHaWAX1q1HbvAqcHrSK2gcrblqIFT/tIyL1RUSAfsB3CZ43I26+edv5JOrXt+VpEZmEe+5cqzSJTMKdxIAB1oHvm2+sH8bNN1s9zJQpST2Fc66SSHS4j8+ALsBU4HFgNtBTVSckuP8m4FzgbexB/2Iw+95QERkabPYG8CMwE/gnwZzfqvo58DI23/e3QZpHJXR1GRJrPolRo2x5WgwfnrZJuBs3trkyxoyBn3+2UWzvvBM2b076qZxzGSRazkkXRKQzNj7UwOQmKTm6d++uEydOzHQyMqNGjdiTaYikdNCnxYth6FDrj9GrFzz5pE0z65yrGkRkkqp2j7Wu1JxFUPRzo4j8n4jcLSKNRKStiPwH+ARYnIoEuwpq0iT28hRPwr3TTvDKK/D001Y81aWL5aiSPQmUcy794hVDPQgcBUzDKppfAT7AiqPyVPWc1CbPldlvv8HGjZa7CItMwn3JJSmtiRaBQYOs7mKffWySpSOOKN6T3TlX9cQLFocBh6rq5cAArHL5FFW9SlWXpjx1ruxuuskGn7rhhm0n4T73XJu8+8ADrdI7hVq2tPk37r8fxo2zecJfeCGlp3TOpVCpdRYislJVG5X0uTLLyjqLH36w7taDBtnQsbG8/DKcdZblPJ58Eo4pdUivpJgxw1pMff45nHiiDaLYtGnKT+ucK6Ny11lgs+L1FZGDROSg4GBbP0eWuUri0kuhbt3S2+gef7z1qmvXDv7wBxuLPMUdJAoKbDrYm26yOo1OnWzwROdc1REvZzGHbXtah6mqtk12opIh63IW771nY3DcdpvNoBTP+vU2Z+v990OPHlZG1KZNypP51VeWy5gyxbp+3HWXzQnunMu80nIW5W46W9llVbDYtMnG3Vi9GqZNs9xFov79bzgzmNvqiSfg2GNTk8aQdetsUMIRIyw+PfWUNbV1zmVWRYqhXFXwz3/aV/URI8oWKMDGHf/qKysrOu44uPDClBdL1a1r08h+8IF9PvBAm3N83bqUntY5VwEeLKq6FSvg6quhd+/y5wratLFKhQsvtLHV998ffvwxmamM6YAD4OuvrTjqzjuhe3eLW865yseDRVV3ww2wfLmN7FeRibJr14aRI6379cyZsNdeVkSVYg0bwsMP2xwfy5fbqLY33WQla865ysODRVU2fbpVUP/5z1BYmJxj/uEP9vV+993hj3+0SSvWr0/OsUtx+OFWknb88ZZR6tXLLs85Vzl4sAhL4bDeKXHppVCvHtx4Y3KPm5cHH35ozWofeMCKpWbNSu45YthhB/jXv+D5563LSNeuFgtTOJyVcy5BHiwi0jSsd9K88w689pp9Dd955+Qfv3Zt6+09ZowFir32sg59aXDiifDtt9CnD5x/Phx6aNGc5s65zPCmsxF5ebGHwGjdGubMSVaykmPTJhulb/16mDoV6tRJ7fnmzrUn+OefwznnlK/VVTmowqOPWgYnJ8dyGYMGVaxqxjlXMm86m4iqNBfqI49Yf4oRI1IfKMAC5vjxNgjhgw/a1HgzZ6b8tCJw9tk2gm3nznDGGVaNksEZc53LWh4sIirFXKgJWL7cerQddFBaxnXaqnZtC06vvmo5rb32StvIgG3b2mCEd94Jr79uw1+NGZOWUzvnAh4sIjI+F2qCbrjB+laMHJmZ8pijjoLJk+2JfdJJ8Le/paU3XU6O1edPmgQtWliXksGDbUR251zqebCICM+FGnHbbWmcCzUB339vxUBnn23lMpnSqpUVS112GTz0EOy7rzVfSoOOHeGzz+Cqq+DZZ21Qwv/9Ly2ndi6rebAIO/VUK2KZPdu+yqahuWiZXHwxNGiQ/Kay5VGrlo3Z8dprVq+z117W5jUNate2W/Dxx9ZyuF8/uOCCbacdd84ljweLWPLy4LTTLKexuJLMHPvmm/a65hrIzc10aooccYQVS3XuDCefbJNwr12bllPvvbf1HzzvPLjvPotXX3xR9brLOFcVeNPZkkyfDnvuacN933pr8hJWHhs32sN482br5ly7dmbTE8vGjdbn4/bbLa0vvmi9wNPk/ffhT3+C+fOhZk1LTkT9+hb3K1OJonOVkTedLY/dd4eBA62O4NdfM5uWhx6y+oq77qqcgQKsWOq226y50s8/26iAzz2XttP362cd+erXLx4owIqnhg9PW1Kcq5Y8WJTmyivh999tyItMWbYMrrvOJjY68sjMpSNRAwZYsVSXLvZVfsiQtBVLNW5ccr3F3Lk21uInn3jdhnPl4cGiNF262AP6nntg1arMpOG666x96N13V52uyy1aWMeIK66wuTb23ttyRmlQUreYnBzrCb7//tCokY27ePbZVjw1efK2uRHnXHEeLOIZPtw6wj3ySPrPPXWqFUENHWptRquSmjXhllusUn7hQiuWevbZlJ+2pO4yTz0FCxbAf/8Lw4bZcFqvvAJ/+YsNWNiokXVMv+ACS+aMGT6AoXNhXsGdiH79bHiN2bPTMiYSYAMj9e9vzXt++AF23DE9502Fn3+2llIffghnnWVNl6Kf6Ek0erTF+J9+spzGzTfHrtxWtTmeJkyw2zxhAnz5ZVExVePGFuN69LBXz57QvHnVyeA5V1Y+B3dF/e9/FjAefNB6LKfD669bEdjIkTaDXVW3aZMVqd1yi/X+fvFFa21WyWzaBN99VxQ8JkywsakikzE1a1YUPCKvpk0zm2bnksWDRUWpWmH3ggX2Lb9WreQctyQbNljXZBFr4pPq86XTO+9YH5bVq62I7fTTM52iuNats+lfwwFk+nT7swAbuyocPPbay2YAdK6qKS1Y1Ex3YqokESvXOPJIK+MYPDi15/vHP6zQ/LXXqlegAJucYvJkOOUUG0Z23DhrbZbCYqmKqlvX6uj33rto2cqVNk5VJIB8+mnRuIo1akD79sUDSOfOlbfVs3OJ8JxFolTtK+OaNVZ/kZOTvGOHLV0K+fn2ZHrzzepbQL5pE1x/vVUo7LknvPSSPWGrsEWLYOLE4jmQpUttXe3a1gIrXP+x++4WWJyrLLwYKlleeglOOMG+Qp5wQnKPHXHOOdby6ptvqvzDMyHvvmu1z6tXW51QqnNtaRSZcDFcgT5pUlEr7O22g27digeQVq2q7/cDV/lVimAhIv2Be4Ec4FFVvS1qvQTrBwBrgMGq+mWwbnvgUaAjoMCZqvppaedLSbDYvNkqZ+vUsaKUZP9XT5lifTv+9jebFi5bLFxoxVLjxlnR1IMP2oCJ1dDmzVbfEQ4gX39t1VRgw35FV6DvtFNm0+yyR8aDhYjkADOAQ4D5wATgZFWdFtpmAHAeFiz2Bu5V1b2DdU8BH6rqoyJSG6ivqitKO2dKggVYg/3Bg20SoKOOSt5xVa08f9Ikq0TPtiY2mzfbULI33AB77GG5uA4dMp2qtFi/3toxhAPItGlFFeitWxcPHt26Wb8Q55KtMgSLfYHrVPWw4PMVAKp6a2ibR4Bxqvqv4PN0oA+wGvgaaKtlSGzKgsXGjVBQYL26Pv00ebmL//s/OPpo64Nw3nnJOWZV9P77Viy1cmVRsVQWlsusWmV9PsIBZPZsWydi8TQcQLp0SV8XIFd9VYbWUM2BeaHP87HcQ7xtmgObgCXAEyLSBZgEXKCqq6NPIiJDgCEArVI1HWqtWjYS7V//WtT/oqI2bLD5rffc03prZ7N+/ayI79RT4cwzYexYax2WZW1RGzaEAw+0V8TSpVaBHgkgb78NTz9t62rVshZX4QDSvn3q2mG47JOuthixvhpG5xJK2qYmsBfwkKp2xXIaw2KdRFVHqWp3Ve2em8o5HwYPhl12Sd6Uq/ffb0VPd99d/ZrKlkezZtYf47rrbOyNHj2snCbL7bijdeq/+mrLiC5caL3UX3nF5sVq3NgG+j3rLAscjRtbsLnkEpuXatasoqIt58oqXTmL+UDL0OcWwIIEt1Fgvqp+Hix/mRKCRdrUrWsTQl9yiQ1jut9+5T/WkiVWTj9ggD0JnMnJgWuvhQMOsMrvnj2tP8aZZ2ZlsVQsItCypb2OO86Wbdli3zvCxVcPPmj1IgA77LBtBfouu2TuGlzVka46i5pYBXc/4GesgvsUVZ0a2uYI4FyKKrjvU9WewboPgT+r6nQRuQ5ooKqXlXbOlNVZRKxebTWPe+9tQ3OU19Ch8Nhj9s15jz2Sl77q5JdfrNd3pD7j4YezrliqIjZutIZ24QAydaq1KQAbJDgcPLp3h+23z2iSXYZkvII7SMQA4B6s6ezjqnqziAwFUNWHg6azDwD9saazf1LVicG+hVjT2drAj8G6UmckSnmwALjpJisT+PJLG7q0rL75xvY77zwbBt2VbPNmG1fquuus0+KLL1pZiyuXNWtsStpwAJk5s2h9QUHxANK1q8137qq3ShEs0i0twWLFCstdHHqoNfUsC1WrzP36a/svbdIkJUmsdsaNsxFsV6yAe++1SSm8WCopfv21eAX6hAk2HBpYqWCnTsUDSMeONhK9qz48WKTS8OE2R/fUqWUbRXXMGDj2WCuHP+eclCWvWlq0CAYNst7fJ59sPd632y7TqaqWFiwoGrrkiy8smERmGa5Xz3Ic4QCy224+hElV5sEilZYsgbw8OP5467CXiPXrrcNZ3brWTNS/npXdli0WpK+5Btq1s5xdly6ZTlW1p2qtqsIB5Msvi2bO3X774nOA9Ojhc4BUJR4sUu2ii6z564wZNl51PHfeCX//uzUPPeSQ1KevOvvgA8tdLF9uxVJDhviTKc02bbIe5+EA8u23xecA6dmzeADZYYfMptnF5sEi1X7+2YLEn/5kLXVKs2iRVdD27m2N5V3FLV5sxVLvvAMnnggHH2yND+JNledSZu1aq44L139Mn160vl27becAqabDgVUpHizSYehQeOIJm6ezefOStzv7bCuumjLFmpy45NiyBW6/Ha680nIW4b/r+vVh1CgPGBn222829Fk4gMwLxmyoUcNKZsMBpFMnnwMk3TxYpMPs2ZZjOP9864kdy+TJ9hXqoovgrrvSl7Zs0qyZ5d6iNW1q3Ztzc4tedeqkP32umEWLihdfTZgAy5bZujp1is8B0qOHzwGSah4s0uX0023shTlz7GEUpgp9+1qrqR9+8F5PqVKjRuJjWjRqVBQ4dtqpeCAJf4689+CScqr27xMOIJMmWR9YsEZv0RXoPgdI8niwSJfvvrO89BVXbDtu1L//DX/8o807ne2DBaZSXp7NOBRt111t0qrFi60F25Ilxd9HPi9dWlQzG2277UoOJLE++zCwSbF5M3z/ffEA8vXX1jMd7FZHV6Cncmi46syDRTodf7y1/587tyj3sG6dDQHaoIF1m/WmsqkzerS1iFqzpmhZWeosVK3DX6xAUlKQiRdcEgksubneRboM1q+3ARDCAeS774oylXl5284B4l1x4qsMQ5Rnj+HDrSjqwQftPdhQHrNnw3vveaBItUhAGD68fK2hRKw3fZMmVkAeTyS4xAss8+ZZh4TFi0sOLg0bll4MFv05i4NLnTpFgSDi99+L5gCJBJDIwAoi1mc2eg4QL1lMnOcsUuGII+Dzzy138fvvVvHdr5/12nbZTdWaBZVUDBbrfaS8JVqDBonVtURe9eun91orgSVLioYwiQSQxYttXa1aFjDCAWTPPbN7DhAvhkq3Tz6B/fe3b6eRsRFGjLAhzZ0rC1WbNTBeXUv4c2RC72gNGpStQr8aBhdVy+RFgseECRZMVq609Q0aWIPFcB1ImzbZU4HuxVDpNnu2tcqJBAqwYSmaNfO2/q5sRGwWo8aNLYcaTyS4xAssCxZYU+7Sgkv9+mWr0K8CvepErGSyVStrbwLWRWfGjOK5jwceKJoDpGnTbecAadYsc9eQKZ6zSIWSWuS0bm3tAp2rLFStqLQsxWKRp2i0SHBJNPdSiYPLhg1Fc4BEAsjUqRZYwOYACec+une3eF7VeTFUupXU1l+k6K/NuapIFVatSrxYbPHikoNLvXplq9Bv0CCj5UGrVxfNARIJILNmFa0vKCgeQAoLq14bBA8W6eY5C+dMJLgkGliWLLGm5rHUrVu2Cv2GDVMeXJYv37YCfeFCW1ezps35EQ4gHTpU7gaRHizSraJt/Z3LVqr2FT5e35bI58WLSw8uZanQT1Jw+fnn4hXoEyZY62oomgMkHEB2263yVKB7sMiE0aPL39bfOZeYcHBJpAPl4sVFk29Eq1OnbBX6222X0FNe1SbDDAePWHOAhANIaWORppIHC+eciwgHl0RyL+ESgrBIcEm03iUUXDZtsgrzcAD55hsb2gRgl122rUBPxxwgHiycc668onMu8epeIqMeRqtdu9TcyvpGucxYkcuX83bi4xm5jJ/ciOkzinIukTlAIkGka9fiDco++tto8kYNZ9fNP7EgpxVzhtxMr3+UrTTDg4VzzqXLmjVlKxYrJbhs2TGX1fVyWSq5zFufy4xfd+LHVbksIZelshMN2+TSomsuBQvGccKnF9KAolzQaurz1V9HlSlgeLBwzrnKau3asvVzWbUq4UPPz2lNi01zEt7ee3A751xlVa9eUbfyRESCS/DSxUvgjNOJVdW+6+afkpZMDxbOOVeVRAUXAeafeTUtNm/bt2tBTitaJOm0PkGhc85VcXOG3Mxqig/8uJr6zBlyc9LO4cHCOeequF7/OJWv/jqK+Tmt2YIwP6d1mSu34/EKbuecc0DpFdyes3DOOReXBwvnnHNxebBwzjkXlwcL55xzcXmwcM45F1e1bQ0lIkuAucCOwNIMJyeTsvn6/dqzVzZff0WuvbWq5sZaUW2DRYSITCypKVg2yObr92vPzmuH7L7+VF27F0M555yLy4OFc865uLIhWIzKdAIyLJuv3689e2Xz9afk2qt9nYVzzrmKy4achXPOuQqq1sFCRPqLyHQRmSkiwzKdnmQTkcdFZLGITAkt20FE3hWRH4KfTULrrgjuxXQROSwzqU4OEWkpImNF5DsRmSoiFwTLs+X664rIFyLydXD91wfLs+L6AUQkR0S+EpHXgs9Zce0iMkdEvhWRySIyMViW+mtX1Wr5AnKAWUBboDbwNdA+0+lK8jUeCOwFTAktuwMYFrwfBtwevG8f3IM6QJvg3uRk+hoqcO27AHsF77cDZgTXmC3XL0DD4H0t4HNgn2y5/uCaLgaeA14LPmfFtQNzgB2jlqX82qtzzqInMFNVf1TVDcDzwDEZTlNSqep4YHnU4mOAp4L3TwF/CC1/XlXXq+psYCZ2j6okVV2oql8G738HvgOakz3Xr6oamYy5VvBSsuT6RaQFcATwaGhxVlx7CVJ+7dU5WDQH5oU+zw+WVXc7q+pCsAcqsFOwvNreDxHJA7pi366z5vqDYpjJwGLgXVXNpuu/B/g7sCW0LFuuXYF3RGSSiAwJlqX82qvzHNyx5i/P5qZf1fJ+iEhD4BXgQlVdKRLrMm3TGMuq9PWr6magUES2B/4jIh1L2bzaXL+IHAksVtVJItInkV1iLKuS1x7YX1UXiMhOwLsi8n0p2ybt2qtzzmI+0DL0uQWwIENpSadFIrILQPBzcbC82t0PEamFBYrRqvrvYHHWXH+Eqq4AxgH9yY7r3x84WkTmYMXLB4nIs2THtaOqC4Kfi4H/YMVKKb/26hwsJgD5ItJGRGoDJwGvZjhN6fAqcEbw/gzgv6HlJ4lIHRFpA+QDX2QgfUkhloV4DPhOVe8OrcqW688NchSISD3gYOB7suD6VfUKVW2hqnnY//X/VPU0suDaRaSBiGwXeQ8cCkwhHdee6Zr9FLcaGIC1kpkFDM90elJwff8CFgIbsW8QZwFNgfeBH4KfO4S2Hx7ci+nA4ZlOfwWvvReWnf4GmBy8BmTR9XcGvgqufwpwTbA8K64/dE19KGoNVe2vHWvd+XXwmhp5rqXj2r0Ht3POubiqczGUc865JPFg4ZxzLi4PFs455+LyYOGccy4uDxbOOefi8mDhsp6IPCkiN2Xo3CIiT4jIryJSpvbvweijB6cqbc6FebBwlU7wEFwUdDqKLPuziIzLYLJSpRdwCNBCVdM+uJ2IXBf0fnauVB4sXGVVE7gg04koKxHJKeMurYE5qro6FelJNRGpzuPLuRAPFq6yuhO4NDKkRZiI5ImIhh9UIjJORP4cvB8sIh+LyEgRWSEiP4rIfsHyeWITRp0Rddgdg0ljfheRD0SkdejYewTrlgcTyJwQWvekiDwkIm+IyGqgb4z07ioirwb7zxSRs4PlZ2FDbO8rIqskmMAoxv5ni03y9LuITBORvWJsU6woTUT6iMj80OfLReTn4BjTRaSfiPQHrgRODM7/dbBtYxF5TEQWBvvcFAmCUfd2OXCdiOwW3LPfRGSpiLwQ6zpc1ebfClxlNREbHO9S4Kpy7L839iBuClyPDTj3f8BuQG/gFRF5RYvmhDgVmx/hc2wimdFAr6Ao7F3gGuBwbJiNd0RkqqpODfY9BRtq5Ehsoq1o/8KGZtgV2AMbKfRHVX1MRDYDf1bVXrEuQkQGAtdh8xNMBNphw7skTER2B84FeqiNVpqHTYAzS0RuAXZTG1sp4ilgEXavGgCvYcNcPxKs3xu7nzth82g8DryDBcraQPeypM9VDZ6zcJXZNcB5IpJbjn1nq+oTasN4v4CNvHmD2iQw7wAbsIdhxOuqOl5V12Nj6ewrIi2xADAnONYmtQmXXgGOD+37X1X9WFW3qOq6cCKCY/QCLlfVdao6GQtigxK8jj8Dd6jqBDUzVXVuGe/FZmymtPYiUktV56jqrFgbisjOWFC8UFVXq41sOhIbsC9igareH9yPtVjwag3sGlzjR2VMn6sCPFi4SktVp2Dfasszf/qi0Pu1wfGilzUMfd46QUyQ21iO5QRaA3sHxVkrRGQFlgtpFmvfGHYFlqvN5hcxl8QnoGmJDQJXbqo6E7gQy6EsFpHnRWTXEjZvjeUWFoau9xGKJtOBba/379i8CV+IzQd+ZkXS6yonDxausrsWOJviD9dIZXD90LLww7s8to75Lzah0g7YuP/zgA9UdfvQq6Gq/jW0b2mjcS4AdogMKx1oBfycYLrmYUVP8aymlPuhqs8FRV2tsfTeHlkV43zrsTmeI9fbSFU7hA8XdexfVPVsVd0V+AvwDxEJ59pcNeDBwlVqwbfiF4DzQ8uWYA/b08SmFj2TxB6opRkgIr3E5j65EfhcVedhOZsCERkkIrWCVw8R2TPB9M8DPgFuFZG6ItIZG0p+dILpehSr6O8W9MnYLVz5HjI5uIYdRKQZlpMArM5CRA4SkTrAOixXtTlYvQjIE5EaQXoXYvUPd4lIIxGpISLtRKR3SQkUkYFic2ID/IoFk80lbe+qJg8Wriq4AatoDTsbuAxYBnTAHsgV8RyWi1kOdMOKmgiKjw7FyuwXAL9g38rrlOHYJwN5wf7/Aa5V1XcT2VFVXwJuDtL3OzAGy/VEewab42AO9rAPt0iqA9wGLA3SvxPWCgrgpeDnMhH5Mnh/OlZRPQ17+L8M7FJKMnsAn4vIKmyynQtUdXYi1+eqDp/PwjnnXFyes3DOOReXBwvnnHNxebBwzjkXlwcL55xzcXmwcM45F5cHC+ecc3F5sHDOOReXBwvnnHNxebBwzjkX1/8DgwPokwbCJX8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = mean_rejection_data\n",
    "x = y.index.values\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (8, 5)\n",
    "plt.xlabel(\"Number of clusters\", fontsize=12)\n",
    "plt.ylabel(\"Rejection rate\", fontsize=12)\n",
    "plt.plot(x, y[\"uniform\"], label=\"Uniform Bootstrap\", color=\"blue\", marker=\"o\")\n",
    "plt.plot(x, y[\"cluster\"], label=\"Cluster Bootstrap\", color=\"red\", marker=\"o\")\n",
    "plt.legend()\n",
    "plt.suptitle(\"Comparison of Rejection Rates\", fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that when the number of clusters is low, it is particularly important to use the cluster robust bootstrap, since rejection with the regular bootstrap is excessive. For a large number of clusters, clustering naturally becomes less important. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}