Chapter11Pscore.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 11: The Central Role of the Propensity Score in Observational Studies for Causal Effects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pandas           : 2.0.3\n",
      "matplotlib_inline: 0.1.6\n",
      "seaborn          : 0.13.0\n",
      "statsmodels      : 0.14.0\n",
      "matplotlib       : 3.7.2\n",
      "numpy            : 1.24.3\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from joblib import Parallel, delayed\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from causalinference import CausalModel\n",
    "font = {'family' : 'IBM Plex Sans Condensed',\n",
    "               'weight' : 'normal',\n",
    "               'size'   : 10}\n",
    "plt.rc('font', **font)\n",
    "plt.rcParams['figure.figsize'] = (6, 6)\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "np.random.seed(42)\n",
    "%load_ext autoreload\n",
    "%autoreload 1\n",
    "\n",
    "%load_ext watermark\n",
    "%watermark --iversions\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[CausalInference](https://github.com/apoorvalal/CausalinferencePy) library for Estimators discussed in Imbens and Rubin (2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhanes_bmi = pd.read_csv(\"nhanes_bmi.csv\")\n",
    "nhanes_bmi.head()\n",
    "z, y, X = (\n",
    "    nhanes_bmi.School_meal.values,\n",
    "    nhanes_bmi.BMI.values,\n",
    "    nhanes_bmi.iloc[:, 3:].values,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regEst(y, z, X):\n",
    "    m0 = sm.OLS(y, sm.add_constant(z)).fit(vcov=\"HC2\")\n",
    "    m1 = sm.OLS(y, np.c_[sm.add_constant(z), X]).fit(vcov=\"HC2\")\n",
    "    m2 = sm.OLS(\n",
    "        y,\n",
    "        np.c_[\n",
    "            sm.add_constant(z),\n",
    "            X - X.mean(axis=0),\n",
    "            z.reshape(-1, 1) * (X - X.mean(axis=0)),\n",
    "        ],\n",
    "    ).fit(vcov=\"HC2\")\n",
    "    res = np.c_[\n",
    "        np.r_[m0.params[1], m0.bse[1]],\n",
    "        np.r_[m1.params[1], m1.bse[1]],\n",
    "        np.r_[m2.params[1], m2.bse[1]],\n",
    "    ]\n",
    "    return pd.DataFrame(res, index=[\"est\", \"se\"], columns=[\"naive\", \"fisher\", \"lin\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>naive</th>\n",
       "      <th>fisher</th>\n",
       "      <th>lin</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>est</th>\n",
       "      <td>0.533904</td>\n",
       "      <td>0.061248</td>\n",
       "      <td>-0.016954</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>se</th>\n",
       "      <td>0.225701</td>\n",
       "      <td>0.220883</td>\n",
       "      <td>0.223635</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        naive    fisher       lin\n",
       "est  0.533904  0.061248 -0.016954\n",
       "se   0.225701  0.220883  0.223635"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regEst(y, z, X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## propensity score stratification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.542402\n",
      "         Iterations 6\n"
     ]
    }
   ],
   "source": [
    "psmod = sm.Logit(z, X).fit()\n",
    "pscores = psmod.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAACUsAAAMMCAYAAABgtmQmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAB7CAAAewgFu0HU+AABrCUlEQVR4nOzde3RddZ3//9dJ0iRN0hYKhXKxKGDhJ+BtoAqilC8WRWFELIOwVBaMBr+CjICMt2GUARERAecr6tQbjBf0WxhFuQl8gSqIQAGlBaUKlMLQG70kzb1N8/uDaaS2hTY57UmyH4+1XCs9e5/PfjerlXyaZ/Yu9fX19QUAAAAAAAAAAGCEq6r0AAAAAAAAAAAAANuCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKYCCu+qqq1IqlVIqlTJ//vwBrdHb25upU6dm9OjROemkk8o7IAAAMCLYewAAAINlXwFAOYilAIa4n//85/1f+G/O/yrhueeey6xZs9LV1ZX/+3//b0VmAAAAkgsvvHCDPcKWfgNhzpw5Of3007PPPvtkzJgxGTNmTPbZZ5+cccYZefTRR7fO4JvJ3gMAALae1atX5z//8z8zderUjB8/PqNGjcr48eNz6KGH5uKLL86KFSs2ey37CgCGMrEUwBDX0tJS6RFe1i677JK3vvWtqaury/Tp0ys9DgAAFNK//Mu/5LzzzhvUGldeeWX+7u/+Lt/4xjcyb968tLW1pa2tLfPmzcuVV16ZN77xjfnWt75Vpom3nL0HAABsHX/5y1+y//775+STT86sWbOyYsWKrFmzJitWrMg999yTz3zmM9lvv/3y+9///mXXsq8AYKirqfQAALy0lStXJkn22GOPXHHFFRWdZVNqamry61//utJjAABAYZ1zzjm57LLLkiTHHntsGhsb86Mf/WiL1rjpppvy8Y9/PH19fdl+++1zyimn5PWvf31KpVIefvjhfP/738+KFSty+umnZ88998yRRx65NX4rL8neAwAAym/BggU55JBDsnTp0lRVVeW4447L4Ycfnu233z5PPfVUrr766sybNy8LFy7MMccck3nz5mX06NEbXcu+AoDhQCwFMMStu7PU7rvvnmOPPbaywwAAAENKX19fzjjjjHzjG99Ikrz//e/PD37wg3z0ox/d4nXOPvvs9PX1ZdKkSfntb3+b3Xbbrf/4Bz7wgZx11ll5y1vekgULFuTss8/OnDlzKvYocAAAoHw++9nPZunSpRk9enSuv/76TJs2bb3jn/rUp/IP//AP+a//+q88++yzueaaa3LqqadusI59BQDDhcfwAQxx62Kp7bffvsKTAAAAQ8natWvzkY98pD+UOuWUU/KjH/0oNTVb/rNxv/rVr/L4448nSb7+9a+v9w2NdXbfffd8/etfT5I8+uijuf322wcxPQAAMFRceeWVed/73pcvfelLG4RSSVJdXZ3Pf/7z/b/e1F2Z7CsAGC7EUgBD3LrH8ImlAACAF+vp6cncuXOTJKeffnq++93vpqpqYP/Uc/PNNydJdt555xx99NGbPO/oo4/OTjvtlCS55ZZbBnQtAABgaBk3blyuvfbanHnmmZs8Z9999+3/eOnSpRs9x74CgOFCLAUwxK27s9R22223Ta7X09OTiy++OAcccEAaGxuz44475ogjjsjMmTNf8n1Tp05NqVTK1KlTNzh21113pVQqpVQq5a677kqS/OAHP8iUKVOy3XbbZcyYMTnssMNyww03vOQ1Ojs785WvfCVvfvObM27cuNTV1WXSpEk59thj88tf/nKgv2UAABiW6uvr88tf/jKXXHJJvv71rw/q0RWzZs1KkhxyyCEvuU6pVMohhxySZNM/Tb657D0AAGBoebm9wDqNjY0bPce+Yn32FQBD15bflx2AbWpbP4bvuOOOy4033tj/646Ojtxxxx2544478oEPfCBXX331gH9afZ0LL7ww55133nqv/frXv86vf/3rfP3rX8/pp5++wXuWL1+eww8/PI888sh6rz/zzDN55plncv311+ekk07KVVddlVGjRg1qPgAAGC4mTJiQc889d9Dr/OUvf0mS7LXXXi977p577pkk+fOf/zyoa9p7AADA8LHu8XpJ8qY3vWmj59hX/JV9BcDQJpYCGOLWPYbvxXeWmj9/fp544olUV1dn3333zcSJE8tyrRtvvDG33357mpub83d/93dpamrK3Llz853vfCdLly7ND3/4w+y55545//zzB3yNZ599NhdddFE+/OEP5+CDD051dXVuvvnm/PSnP02SnHPOOTn66KOzxx57rPe+z3zmM/2bioMOOigf+MAHstNOO+W///u/c8011+TBBx/Mj3/844wZMybf+ta3Bv5JAACAglm6dGna29uTJE1NTS97/pgxY5K88IMdy5cvz/jx47f4mvYeAAAwvHzpS19K8sIj+z70oQ9tcNy+wr4CYDgRSwEMcS++s9TChQtz8skn57bbbus/XiqVcvjhh+eyyy7L6173ukFd6/zzz8/cuXOz9957r/f6Oeeck8MPPzxz5szJJZdckubm5uy2224DvsasWbNy0EEH9b928sknZ//99895552X7u7u/PCHP8znPve5/uNr167t33hMmTIlv/3tb1NdXd1//Oyzz84XvvCFXHnllTn55JMHNBcAABTV888/3//xi7/O3pSamr/+c9KSJUsG9E0New8AABj6Ojo68oc//CGXXXZZrr322owdOzYzZ87MhAkTNjjXvsK+AmA4Gdw9BwHY6tbdWaqxsTHTpk1bL5RKkr6+vtxxxx15y1veknvuuWdQ17riiis22FQkyQ477JDvfe97SZKurq5cc801A77GCSecsN6mYp1//ud/7n/U4O9+97v1jj3//PP90dj06dM32GiVSqWcf/75eeyxx3LwwQcPeDYAACiizs7Obf5eew8AABiarrjiipRKpZRKpTQ2NuaQQw7JzTffnBNOOCH33Xdfpk2bttH32Vf8lX0FwNAnlgIY4tZ9QX355Zdn1113zW9+85u0tLSkra0t9913X6ZPn54kaW9vz0knnTSoDcmUKVM2eezAAw/MvvvumyS55ZZbBnyNQw45ZKOv19bW9l9/wYIFm3z/Sz27e6eddhrwXAAAUFSD2UN0dHQM6H32HgAAMHyMGzcuu+22W0ql0ibPsa/YOPsKgKFJLAUwhK1ZsyaTJ09OU1NTdt9999x888059NBDM3bs2DQ2NmbKlCmZOXNmTjrppCQvfEH+k5/8ZMDXe/FtbzfmNa95TZL0P2d7IBoaGjZ5bN2te9va2jZ4fd2x6667Lr29vQO+PgAAsL6+vr4Bv/elvlnyUuw9AABgaDrmmGPys5/9LD/72c9yzTXX5Pzzz09DQ0Muu+yyHHDAAfnOd76z0ffZV6z/un0FwNAmlgIYwmpqajJnzpysWrUqP/nJTzb5nO8LLrig/+Mbb7xxq80zZsyYJMnSpUvT3d1d9vXX/f7+duNQKpXysY99LEly99135w1veEMuv/zyzJ49e8A/cQIAALzgpf7x/+WMHj26jJP8lb0HAABUxl577ZVjjz02xx57bN7//vfnX//1X/Poo4/muOOOy+rVq/PRj350g8fOJfYVL2ZfATD0iaUAhomqqk3/X/aee+6ZPfbYI0ny0EMPbZN5Wltbt8l11jnvvPPyT//0T6murs6cOXNy9tln56CDDsrYsWNz0EEH5Stf+coGP70BAAC8vMF8Y2JrfVPjxew9AACgsmpra/Ptb387o0ePTm9vb7785S9vcI59xfrsKwCGNrEUwAix6667JkmWLVu2Ta430NviDlR1dXWuuOKK/PnPf86FF16YadOmZYcddkhvb29mz56df/7nf87++++fv/zlL9t0LgAAGO523HHH/o835/EQa9as6f94p5122iozvZi9BwAAVN748ePz5je/OUny//7f/9vgsXv2FeuzrwAY2sRSACNEbW1tkvU3GFvTutvXbmuvetWr8rnPfS633nprnn/++Tz++OP57Gc/m+rq6jz99NM55ZRTKjIXAAAMVxMmTEhjY2OSbNZPNq9atSpJMnbs2IwfP36rzpbYewAAwFCx7oe2V61alZUrV653zL5i4+wrAIYmsRTACNHe3p4k2X777bfaNdZtXiZMmJC6urqtdp0tMXny5Hzxi1/MJz7xiSQvPP/7v//7vys7FAAADDN77713kuSJJ5542XOfeuqpJC98Lb612HsAAMDQtnbt2g1es694efYVAEODWApgiGttbc1DDz30suet23zss88+A77Wy92V6rHHHkuSvPa1rx3wNbaWqVOn9n/89NNPV24QAAAYhg477LAkyb333rvB4zRerK+vL7/97W+TJG9729sGfD17DwAAGDo6Ozvz5JNPvux5zz33XJJk1KhRG/3BbfuKzWdfAVBZYimAIezRRx/NQQcdlOOOO+4lb1t71113ZcWKFUn+uhkZiAceeGCTxx588MH86U9/SpK84x3vGPA1Burb3/52fvzjH2/y+PPPP9//8bpb/QIAAJvnqKOOSpIsWrQoN9100ybPu+mmm7J48eIkg9sX2HsAAMDQ8OSTT+aQQw7Je97znnR3d2/yvOXLl+d3v/tdkuRNb3pTqqo2/DazfcVf2VcADG1iKYAh7J577sm8efPy9NNP513veldaWlo2OGfFihU588wzkyQ1NTU59dRTB3y9T3ziE/nLX/6ywevLli3rf252fX19TjzxxAFfYyC+9rWvpbm5OaeeemquueaaDY6vXr063/zmN/vnG8zdtQAAoIiOPPLI/sdfnHHGGf0/Mf5izz33XM4444wkyWte85pMmzZtwNez9wAAgKHhpz/9aX7/+99n7ty5mT59ejo7Ozc4p6enJx/96Ef7j637mv1v2Ve8wL4CYOirqfQAAGxac3Nz7rjjjvz0pz/Nb37zm+y555455ZRTsv/++6euri5z587N9773vSxatChJ8vnPfz677777gK933nnnZf/998/JJ5+cAw88ME1NTZkzZ06+853vZOnSpUmST37yk4O6xkCcdNJJueyyy7JgwYKcdNJJ+fa3v533vOc92XnnnfPss8/mu9/9bv9PiJx88smpr6/fpvMBAECl/PnPf86sWbM2eP3xxx/v//inP/1pdthhh/WOH3bYYXn1q1/d/+uqqqp89atfzd///d9n/vz5OeCAA3Lqqafm9a9/fZLk97//fb73ve9l+fLlqaqqyqWXXppSqTTgue09AABgaPjUpz6V++67L9dff31uuOGG7LXXXjnllFMyefLk1NXVZd68efnRj36UefPmJUkOPfTQfOhDH9roWvYVL7CvABj6xFIAQ9yPf/zjTJo0KZdeemmWL1+er371qxucUyqV8slPfjKf+9znBnWto48+OrfccktmzJiRGTNmbHD8xBNPzPnnnz+oawzEhAkT8pvf/CbTp0/PAw88kDvvvDN33nnnBucdeOCB+cpXvrLN5wMAgEq555578pGPfOQlz/n0pz+9wWvf//7314ulkhf2A1/72tdyzjnnZPny5bn00ks3eN+oUaNy+eWX9z9eY6DsPQAAYGioqqrKddddlzPPPDPf/OY3s3Dhwlx00UUbPXfq1KmZOXNmamo2/S1m+wr7CoDhwGP4AIa4qqqqXHLJJbn//vtz6qmnZq+99kpdXV0aGhqy99575x//8R8ze/bsXHLJJYP6CYx1/uu//isXXXRR9ttvvzQ0NGT8+PGZOnVqfvzjH+fHP/7xRp9Dvi1MmjQp99xzT374wx9m2rRpecUrXpFRo0Zlu+22y1ve8pZ87Wtfy913350xY8ZUZD4AABgJPv7xj+f+++9Pc3NzXv3qV6exsTGNjY159atfndNOOy2zZ8/O6aefXpZr2XsAAMDQUF1dnSuvvDKPPvpompubs/fee6e+vj719fV5xStekfe973257rrrcscdd2THHXd82fXsK+wrAIa6Ul9fX1+lhwAAAAAAAAAAANja3FkKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIVQU+kB/lZXV1fmzJmTJJkwYUJqaobciAAAMGStWbMmS5cuTZIccMABqa+vr/BElWFfAQAAA2df8QL7CgAAGLihvK8Ycl/Zz5kzJ1OmTKn0GAAAMOzdf//9Oeiggyo9RkXYVwAAQHnYV9hXAADAYA21fYXH8AEAAAAAAAAAAIUw5O4sNWHChP6P77///uyyyy4VnAYAAIaXhQsX9v/k84u/ti4a+woAABg4+4oX2FeMPDNmzEjbs39K05rn03z04O8aNuOG+9NWs2Oadt83zc3NZZgQAGDkGMr7iiEXS734md+77LJLdt999wpOAwAAw9eLv7YuGvsKAAAoD/uKF9hXjAzbbbddqloaM3Z1e3afMG7w6zWNTtWoxozdbjt/PgAAXsJQ21d4DB8AAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEGoqPQAAAAAAAADASNPX15f29va0tramq6srvb29lR4Jtkh1dXXq6+szduzYNDY2plQqVXokgLIQSwEAAAAAAACU0dq1a7NgwYJ0dnZWehQYsDVr1qS7uzstLS0ZPXp0Jk2alKoqD68Chj+xFAAAAAAAAECZ9PX1bRBKlUqlVFdXV3Aq2HK9vb3p6+tLknR2dmbBggXZY4893GEKGPbEUgAAAAAAAABl0t7e3h9KVVdXZ+LEiWlqanJHHoadtWvXpq2tLYsWLUpvb286OzvT3t6epqamSo8GMCj+iwwAAAAAAABQJq2trf0fT5w4MWPHjhVKMSxVVVVl7NixmThxYv9rq1atquBEAOXhv8oAAAAAAAAAZdLV1ZXkhUfvuQMPI0FTU1P/o/de/HhJgOFKLAUAAAAAAABQJr29vUleeASfO0oxElRVVaW6ujrJX/98Awxn/usMAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAGCEmzp1akqlUq666qpKjwJQUTWVHgAAAAAAAACgqC6/bV6lR9hqzpo2udIjAMAG3FkKAAAAAAAAgGHriiuuyBe+8IVKj7HZfv7zn+cLX/hC5s+fX+lRAApJLAUAAAAAAADAsHXFFVfk/PPPr/QYm+3nP/95zj//fLEUQIWIpQAAAAAAAAAAgEIQSwEAAAAAAAAwLD377LNZsmRJpcfYbJ2dnfnTn/5U6TEACk0sBQAAAAAAAMCwctddd6VUKuUVr3hFOjs7kySlUimlUilTp07tP2/da/Pnz88DDzyQI444IuPGjcuxxx673nqrV6/ON7/5zbzlLW/JuHHjMnr06LzmNa/Jeeedl5aWlk3OcfPNN+cd73hHdt5559TV1eVVr3pVTjvttCxYsGCDc0ulUhoaGnLfffclSQ4//PD++f7WQOd55plncuqpp2bXXXdNfX199txzz5x77rn9nyMAkppKDwAAAAAAAAAAW6Kmpibjxo3L2rVrs2rVqiTJuHHjkiRNTU0bnP/UU0/l2GOPTWtra8aOHZvGxsb+Yx0dHTnmmGNyxx13JHkhaqqurs4f//jHXHjhhZk5c2ZmzZqVnXfeeb01P/e5z+Wiiy5KklRXV6eqqirz58/PjBkzct111+W3v/1tJk+e3H/+uvlaW1vT19eXxsbG1NRs+C37gc7zpz/9KYceemiWLVuWJBk1alQWLFiQSy+9NPPnz9/8Ty7ACOfOUgAAAAAAAAAMK4ceemhWrlyZRx55pP+1lStXZuXKlbnhhhs2OP/CCy/Mv//7v2fVqlVpaWnJj370o/5jn/nMZ3LHHXdk7733zk033ZRVq1alo6MjDzzwQKZMmZLHH388Z5xxxnrr3XnnnbnooosyefLk3HHHHenu7k5XV1cefPDBvP71r8+yZcvyyU9+cr33rJtv0qRJSZIbbrih/7UXG8g8SfKhD30oy5Ytyxve8Ib89re/TXd3dzo6OvKrX/0qjzzySGbNmrXFn2eAkUgsBQAAAAAAAMCI9t73vjcnn3zyBnedWrx4cb75zW+mrq4uN998c4466qg0NjZm1KhROfDAA/PLX/4yjY2Nue666/LUU0/1v+9tb3tb/uVf/iV33nlnDj/88P47S73xjW/Mt771rSTJrbfemtWrV2/RnAOdZ9asWXnggQcyevTo3HDDDTn44INTKpVSW1ubI488Mrfcckvq6+sH8RkEGDnEUgAAAAAAAACMaIcddthGX//FL36R1atX54gjjsjee++9wfGddtopb3/729PX17feHauqq6tzwQUXZNddd93gPa973euSJN3d3Vm6dOkWzTnQeX75y18mSY4//viNzvSqV70qb3jDG7ZoFoCRSiwFAAAAAAAAwIg2ZsyYjb4+e/bsJNlomLTOK1/5yiTJ73//+8261ovv4NTT07N5Aw5ynj/84Q9Jkn322WeT76utrd2iWQBGqppKDwAAAAAAAADDzb2PLUhr79LUz19e1nWbmprS3Nxc1jWBTVu0aFGS5D/+4z9y9dVXb/Scrq6uJNngLlHd3d2ZMWNGrr/++vz5z3/OqlWrsnbt2orMs+5948ePH9T1AYpALAUAAAAAAABbqGf1mpR621PbsSitT88pz6L145JMKs9awGZpb29P8kL41N3d/ZLntrW19X+8atWqHH744XnwwQeHxDwdHR1lnQNgJBNLAQAAAAAAwADUl1Zn3NqVGbt6yaDXWtXRnb4yzARsmcbGxiTJpz71qVx88cWb/b5/+7d/y4MPPpgdd9wx/+f//J8cfvjhmTBhQqqqqpIkpVJpm87T0NAwoOsBFJFYCgAAAAAAAAaooa4mZx9/6KDXuWzm3WktwzzAltl5552TJJ2dnVv0vmuvvTZJ8sUvfjHvf//7Kz7PxIkTM3fu3CxfXt5HgwKMRFWVHgAAAAAAAAAAKuHAAw9Mkvzxj3/c5DkrVqzIokWL1nvt2WefTZJMnjx5g/Off/75bT7Pa1/72iTJ448/vsn39fT0DHgugJFELAUAAAAAAABAIf393/99Ro0alTvvvDPz5s3b6Dmf/vSnc8ABB+Sxxx7rf62pqSlJsnDhwg3O/973vrfN5znmmGOSJDNnztzoTE8//XQefvjhAc8FMJKIpQAAAAAAAAAYlqqqBvct74kTJ+a0007LmjVr8s53vjO/+tWvsmbNmvT19eUvf/lLPvzhD2fGjBnZfvvts9dee/W/b8qUKUmSCy+8MHPnzk2SLF68OBdccEGuueaaAc880HkOO+ywHHjggens7MwxxxyT++67L319fenp6cntt9+eCy64IG9605sG86kCGDHEUgAAAAAAAAAMSzvttFNqamqSvBAaTZw4Meecc84WrXHxxRdn6tSpeeqpp/LOd74zo0ePTn19fV796lfnu9/9bnbbbbf84he/SF1dXf97PvWpT6WqqiqPPfZYDjjggNTV1WXixIn5t3/7t3zta1/LxIkTk7zweLyLL754vevtuuuuSZLjjjsuEydOzMEHHzzoeUqlUn7wgx9k/PjxefDBB/PmN7859fX1aWhoyJFHHpnjjjtuiz4nACOZWAoAAAAAAACAYam+vj4XXXRRxo8fn8WLF2fx4sVpaWnZojUaGxtz66235sorr8zBBx+c0aNHZ9SoUdlvv/3y2c9+Ng8//HD23Xff9d7zv/7X/8pNN92Ut771rWlsbExtbW3e+ta35pZbbsnb3va2XHLJJZkwYUJWrVqVrq6u9d77hS98Ia985SuzcuXKLF68OEuXLh30PEmy77775qGHHsqHPvSh7LzzzqmqqsprX/vaXHvttXnXu961RZ8TgJGsptIDAAAAAAAAABTVWdMmV3qEYe/cc8/Nueeeu9FjfX19m7XGqFGj8rGPfSwf+9jHNvu673jHO/KOd7xjo8c++MEP5oMf/OBGj7397W/PU089VfZ5kmSPPfbI1VdfvdFjd9111xatBTBSubMUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAh1FR6AAAAAAAAAIAimjFjRtra2io9xlbV1NSU5ubmSo9RWK985Svz9NNP584778zUqVMrPQ7AkCCWAgAAAAAAAKiAtra2tC5ekHS1VHqUraN+XJJJlZ4CANYjlgIAAAAAAAColK6WlFqeyZiGukpPUlarOrrTV+khhogvfOEL2W677fKJT3yi0qMAELEUAAAAW0kRHiXwtzxaAAAAgIEY01CXs48/tNJjlNVlM+9Oa6WHGCLOP//87LHHHmIpgCFCLAUAAMBWMeIfJfC3PFoAAAAAAGDIE0sBAACw9YzQRwn8LY8WAAAAADbmoYceqvQIAPwNsRQAAABb1Uh8lMDf8mgBAAAAqJwlS5bky1/+cm644YYsWLAgY8eOzT777JOPf/zjmT59ekql0nrnP/bYY7nkkktyxx13ZPHixRk7dmze9KY35Ywzzsg73/nODda/6667cvjhh2ePPfbI/Pnz87Of/Sxf/OIX88c//jGjR4/OMccck0suuSQTJkzof89VV12VU045pf/XTz/9dP8cJ598cq666qokyfz58/OqV70qSdLX15dbb701//qv/5o5c+bkIx/5SK644or+Ndrb2/Pv//7vmTlzZubNm5ckmTx5cv7hH/4hZ555ZhoaGsry+QQY6cRSAAAAAAAAAAxLjz76aN7+9rdn0aJFSZL6+vosXbo0S5YsyW9+85ucdtpp+eY3v9kfKl1//fU54YQT0t3dnSRpaGjIsmXLcuONN+bGG2/Mpz/96XzpS1/a5PVuv/32TJ8+PVVVVent7U1HR0euuuqqzJ49O7Nnz05d3Qt3166trc24ceOyZs2atLe3p1QqZezYsf3X3Jjf/e53Ofroo7NmzZqMHTs2o0eP7j+2dOnSHHHEEZkzZ06SZNSoUSmVSnn44Yfz8MMP55prrsntt9++XrAFwMZVVXoAAAAAAAAAANhSa9asyfvf//4sWrQo73rXu/KnP/0pnZ2d6ezszLXXXpvx48fnP/7jP/KDH/wgSbJgwYJ88IMfTE9PTz71qU9l8eLFaW9vz6pVq/Ld7343TU1Nufjii/Pzn/98o9fr7OzMt771rfziF79Ie3t7Vq5cmUsvvTRVVVWZO3durr766v5zTzrppKxcuTI33HBDkmTSpElZuXJlVq5cmW984xsbXf/iiy/OzJkz+9d+cbT14Q9/OHPmzMmb3/zm3Hfffenu7k5XV1fuu+++TJkyJY888kg+8pGPlOkzCzCyiaUAAAAAAAAAGHZ+9atfZe7cuXnFK16R6667Lvvss0+SpK6uLu973/vywx/+MGPHjk1ra2uS5NJLL82qVavykY98JBdffHF22mmnJEljY2NOPfXU/ojpggsu2Oj1lixZkjPPPDPvfve7U1tbm7Fjx+acc87JcccdlyS5+eabB/X7+fCHP5z3vOc9691RKkkefvjh/OIXv8iECRNyyy23ZMqUKSmVSimVSpkyZUpuueWW7Ljjjrn++uvzyCOPDGoGgCIQSwEAAAAAAAAw7Nxyyy1Jkve+972pr6/f4PhRRx2VZ555JmeccUaS5LrrrkuS/OM//uNG1zvxxBPT1NSUhx56KE8++eQGx5uamvLWt751g9ePPvroJMkTTzwxsN/I/zjssMM2+vq1116bJJk+fXrGjRu3wfHtt98+06dPX+9cADZNLAUAAAAAAADAsDN37twkyZ577rnJc8aOHZskWbhwYZ577rkk6b8D1d+qqanJXnvtlSR54IEHNji+ww47pFQqbfD6brvtliRpa2vbguk3NGbMmI2+Pnv27CSbnvvFxzY2NwDrE0sBAAAAAAAAMOw8//zzSV547N7LWbRoUf/HL3X+ujtULVy4cLPnqKmpSZKsXbt2s9+zJdbNXu65AYpKLAUAAAAAAADAsNPR0bHZ57a3t2/R2oO9S1Q5bcnsQ2lugKFKLAUAAAAAAADAsNPQ0LDZ5zY2Nm7R2k1NTVs6zlazJbMPpbkBhiqxFAAAAAAAAADDzo477pgk6e7uftlzd9555/6PX+r8rq6uJMnEiRMHOV35rJt9uM0NMFSJpQAAAAAAAAAYdvbff/8kyZNPPrnJc+bNm5ck2XXXXbPLLrskSR5//PGNnrtmzZo88cQTSZKDDjqonKMOyoEHHphk03O/+NhQmhtgqBJLAQAAAAAAADDsvPOd70yS/OxnP+u/s9KLPfHEE3nd616X008/PUnyvve9L0ny3e9+d6Pr/eQnP0lbW1te//rXZ6+99tpKU2+56dOnJ0muu+66tLS0bHB8xYoVufbaa9c7F4BNE0sBAAAAAAAAMOy84x3vyP77759nnnkm06dP77+LVE9PT2666aYcccQR6erqym677ZYk+eQnP5mmpqZ8+9vfzmc+85ksWbIkSdLR0ZHvf//7+d//+38nSc4777yyzVhVNfhvyb/xjW/M0UcfnSVLluSoo47KAw88kL6+vvT19eWBBx7IUUcdleeffz7HHHNMXve616333lKpNOjrA4w0YikAAAAAAAAAhp2ampr85Cc/yc4775wbb7wx++yzT0aPHp36+vq8+93vztNPP50TTjgh5557bpJkjz32yH/+539m1KhRufjii7PzzjunsbExTU1NOfXUU9PW1pZzzz03xx13XNlm3HXXXZMkzzzzTCZOnJiJEyfmiiuu2OJ1vvOd72T//ffPvffemylTpvT/PqdMmZL77rsv+++/f7797W9v8L51jx48+uijs9122230zlQARVNT6QEAAAAAAAAAimxVR3cum3l3pccoq1Ud3cm4rX+d/fbbL4888ki+/OUv55e//GUWLFiQHXfcMfvuu29OO+20nHjiievd3em9731vHnrooVxyySW54447snjx4uywww5505velNNPPz1HHXVUWefbe++980//9E/5/ve/n8WLFydJ2tratnidnXfeOb/73e/yta99LTNnzuy/i9Z+++2X448/PmeeeWYaGxs3eN/nP//5nHbaaVmwYEH/3aiAwbn8tnllXe+saZPLuh4vTywFAAAAAAAAUCn149KXpLXSc5TbuCT126CWSrLTTjvlq1/9ar761a9u1vn77bdfrr766s1ef+rUqS8ZGb3c8SuuuGKTd5N65StfudkBU2NjYz772c/ms5/97Gadn7zwqML58+dv9vkARSCWAgAAAAAAAKiApqamJJMqPcZW9cLvEQCGDrEUAAAAAAAAQAU0NzdXegQAKJyqlz8FAAAAAAAAAABg+BNLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAACUSVXVC9+C7e3tTV9fX4WngcHr6+tLb29vkr/++QYYzvw/GQAAAAAAAECZ1NbWJnkhMOnu7q7wNDB43d3d/eHfuj/fAMOZWAoAAAAAAACgTBobG/s/bm1treAkUB4v/nP84j/fAMOVWAoAAAAAAACgTJqamvo/XrZsWZYtW9b/CDMYTnp7e/v/DK/z4j/fAMNVTaUHAAAAAAAAABgpamtrM2HChCxdujRJsmTJkixZsiTV1dUplUoVng42T19f3waR34QJEzyGDxgRxFIAAAAAAAAAZbTDDjukp6cnLS0t/a+5uxTD2bhx47LDDjtUegyAshBLAQAAAAAAAJRRqVTKrrvumvHjx2flypXp6OgQSzHsVFdXp6GhIdttt13q6+srPQ5A2YilAAAAAAAAALaC+vr6TJw4sdJjAAAvUlXpAQAAAAAAAAAAALYFsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKoabSAwAAAAAAAAAAwNZy+W3zKj0CQ4g7SwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKoabSAwAAAAAAAEDR3fvYgrT2Lk39/OVlW7OpqSnNzc1lWw8AYCQQSwEAAAAAAECF9axek1Jve2o7FqX16TmDX7B+XJJJg18HAGCEEUsBAAAAAADAEFBfWp1xa1dm7Oolg1pnVUd3+so0EwDASCOWAgAAAAAAgCGioa4mZx9/6KDWuGzm3Wkt0zwAACONWAoAAAAAAAAAYCu5/LZ5ZV3vrGmTy7oeFE1VpQcAAAAAAAAAAADYFsRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIVQU+kBAAAAAABgUy6/bV6lR9gmzpo2udIjAAAAFII7SwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAh1FR6AAAAAAAAAAAANs/lt82r9Agbdda0yZUeATaLO0sBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKISaSg8AAAAAAFBOl982r9IjbBNnTZtc6REAAABg2HFnKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCDWVHoChacaMGWlra6v0GGwFTU1NaW5urvQYAAAAAAAAAADbnFiKjWpra0vr4gVJV0ulR6Gc6sclmVTpKQAAAAAAAAAAKkIsxaZ1taTU8kzGNNRVehLKYFVHd/oqPQQAAAAAAAAAQAWJpXhJYxrqcvbxh1Z6DMrgspl3p7XSQwAAAAAAAAAAVFBVpQcAAAAAAAAAAADYFsRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAg1lR4AAAAAAAAANmbGjBlpa2sry1qzZs1Kbcei7F7fWZb1AAAYnsRSAAAAAAAADEltbW1pXbwg6WoZ9Fpdyxem1Nee3lFlGAwAgGFLLAUAAAAAAMDQ1dWSUsszGdNQN8h1Vv7Pd8bUUgAARSaWAgAAAAAAYEgb01CXs48/dFBrzPrDk8na3jJNBADAcFVV6QEAAAAAAAAAAAC2BbEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEKoqfQAAAAAAAAAAABDyeW3zav0CMBW4s5SAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAAqhptIDAAAU0eW3zav0CJTRWdMmV3oEAAAAAEa4GTNmpK2traxrNjU1pbm5uaxrAsBQJ5YCAAAAAAAAGOLa2trSunhB0tVSngXrxyWZVJ61AGAYEUsBAAAAAAAADAddLSm1PJMxDXWDWmZVR3f6yjQSAAw3YikAAAAAAACAYWJMQ13OPv7QQa1x2cy701qmeQBguKmq9AAAAAAAAAAAAADbglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACiEmkoPAAAAAAAAAMC2c+9jC9LauzT185eXdd2mpqY0NzeXdU0AKDexFAAAAAAAAECB9Kxek1Jve2o7FqX16TnlWbR+XJJJ5VkLALYisRQAAAAAAABAwdSXVmfc2pUZu3rJoNda1dGdvjLMBADbglgKAAAAAAAAoIAa6mpy9vGHDnqdy2bendYyzAMA20JVpQcAAAAAAAAAAADYFsRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQaio9AAAAAMBQM2PGjLS1tVV6jG2qqakpzc3NlR4DAAAAALYqsRQAAADA32hra0vr4gVJV0ulR9k26sclmVTpKQAAAABgqxNLAQAAAGxMV0tKLc9kTENdpSfZqlZ1dKev0kMAAAAAwDYilgIAAADYhDENdTn7+EMrPcZWddnMu9Na6SEAAAAAYBupqvQAAAAAAAAAAAAA24JYCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEIQSwEAAAAAAAAAAIUglgIAAAAAAAAAAApBLAUAAAAAAAAAABRCTaUHAAAAAABgy11+27xKjwAAAADDjlgKAIYJ/wgOAAAAAAAAMDhiKQAAAAAAAIAkM2bMSFtbW1nXbGpqSnNzc1nXBAAGTiwFAAAAAAAAkKStrS2tixckXS3lWbB+XJJJ5VkLACgLsRQAAAAAAADAOl0tKbU8kzENdYNaZlVHd/rKNBIAUD5iKQAAAAAAAIAXGdNQl7OPP3RQa1w28+60lmkeAKB8qio9AAAAAAAAAAAAwLYglgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACkEsBQAAAAAAAAAAFIJYCgAAAAAAAAAAKISaSg8AAAAAAAAAAMDwdvlt88q63lnTJpd1PVjHnaUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUglgKAAAAAAAAAAAoBLEUAAAAAAAAAABQCGIpAAAAAAAAAACgEGoqPQAAAAAAAAAAI9flt80r21pnTZtctrUAKCZ3lgIAAAAAAAAAAApBLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKoabSAwDbxr2PLUhr79LUz19e6VHYCpqamtLc3FzpMQAAAAAAAABgSBNLQUH0rF6TUm97ajsWpfXpOZUeh3KqH5dkUqWnAAAAAAAAAIAhTywFBVJfWp1xa1dm7OollR6FMlnV0Z2+Sg8BAAAAAAAAAMOEWAoKpqGuJmcff2ilx6BMLpt5d1orPQQAAAAAAAAADBNiKQAAAAAAAAAG7N7HFqS1d2nq5y/f6PFZjy8d0Lp1oxty8LtPGMxoALABsRQAAAAAAAAAA9azek1Kve2p7ViU1qfnbHC8dsWWPyeju3pMurNLOcYDgPWIpQAAAAAAAAAYlPrS6oxbuzJjVy/Z4NgupfYtWqujqyctdUmPWAqArUAsBQAAAAAAAMB67n1y2Wafu7y9J1XpTVPdqBz8d//fBscP3sJrX3Pr7LT0beGbAGAzVVV6AAAAAAAAAAAAgG1BLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCqKn0AAAAAAAAAAAMzr1PLtvsc5e396QqvWmq69ui98FQd/lt8yo9AjAMuLMUAAAAAAAAAABQCGIpAAAAAAAAAACgEMRSAAAAAAAAAABAIYilAAAAAAAAAACAQhBLAQAAAAAAAAAAhSCWAgAAAAAAAAAACqGm0gMAAAAADEf3Prms0iOUxfxl7VnY15qerqUp3TZvg+NnTZtcgakAAAAAYOtwZykAAAAAAAAAAKAQ3FkKAAAAAAAAAAAq4PKN3Ol7oNwhfPO4sxQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFUFPpAQAAAAAAAABGmnsfW5DW3qWpn7+8LOvNmjUrtR2Lsnt9Z1nWA4CiEksBAAAAAAAAlFnP6jUp9bantmNRWp+eM+j1upYvTKmvPb2jyjAcABSYWAoAAAAAAABgK6gvrc64tSszdvWSwS/WtfJ/vrurlgKAwRBLAQAAAAAAAGwlDXU1Ofv4Qwe9zqw/PJms7S3DRABQbFWVHgAAAAAAAAAAAGBbEEsBAAAAAAAAAACFIJYCAAAAAAAAAAAKQSwFAAAAAAAAAAAUQk2lBwAAAABg6Lr8tnmVHmGbOGva5EqPAAAAAMA24M5SAAAAAAAAAABAIbizFAAAAAAAADAszZgxI21tbWVbb9asWantWJTd6zvLtiZQLPfe+NN0d3aUZa35f/x9elf3pGbUqOzx/72hLGsmSd3ohhz87hPKth4MN2IpAAAAAAAAYFhqa2tL6+IFSVdLWdbrWr4wpb729I4qy3JAAXV3dqSvdWHqelcNfrHW/07t2jVprB+V2hWNg18vSXf1mHRnl7KsBcOVWAoAAAAG6d7HFqS1d2nq5y+v9CjbRFNTU5qbmys9BsAm7yQx6/GlFZhm2/AT4ACwEV0tKbU8kzENdWVYa+X/fAdVLQUMXF3vqozrXpiG+tpBrbNgTWvqa5IJNVXZpbRs0HN1dPWkpS7pEUtRcGIpAAAAGKSe1WtS6m1PbceitD49p9LjbF3145JMqvQUAEk2fSeJ2hWtFZpo6/IT4ACwaWMa6nL28YcOep1Zf3gyWdtbhomAomuor82JRx44qDUefvyZJL2pr60e9FpJcs2ts9PSN+hlYNgTSwEAAEAZ1JdWZ9zalRm7ekmlR9lqVnV0x7+nAUPORu4ksUupvYIDbR1+AhwAAADKQywFAAAAZdJQV1OWn2Qeqi6beXdG5r1agOHub+8kce+Tg388xVDjJ8ABYGTakq9blrf3pCq9aarrG5Ff70Cl3HvjT9Pd2VG29Z545P5s3/t8Gpq6y7YmUF5iKQAAAAAAAACgkLo7O9LXujB1vavKsl5159Kk1JnetWVZDtgKxFIAAAAAAAAAQGHV9a7KuO6FaaivHfRaC9a0pqomSUYNei1g6xBLAQAAAAAAAACF1lBfmxOPPHDQ6zz8+DNJegc/ELDVVFV6AAAAAAAAAAAAgG1BLAUAAAAAAAAAABSCWAoAAAAAAAAAACgEsRQAAAAAAAAAAFAIYikAAAAAAAAAAKAQxFIAAAAAAAAAAEAhiKUAAAAAAAAAAIBCEEsBAAAAAAAAAACFUFPpAQAAAIDh4d7HFqS1d2nq5y+v9Chb3axZs1LbsSi713dWehQAAAAAoIzEUgAAAMBm6Vm9JqXe9tR2LErr03MqPc5W1bV8YUp97ekdVelJAAAAeLHLb5tX1vXOmja5rOsBMPSJpQAAAIDNVl9anXFrV2bs6iWVHmXr6lr5P/9qMvJrqblPPJdla5and/TKSo+yzdSNbsjB7z6h0mMAAAAAUAFiKQAAAGCLNNTV5OzjD630GFvVrD88maztrfQY28TqNb2p7e1Mfe/zqV1R3p/QHoq6q8ekO7tUegwAAAAAKkQsBQAAAFBwDdWrM6FmVXYpLav0KFtVR1dPWuqSHrEUAAAAQGGJpQAAAABIfW11TjzywEqPsVVdc+vstPRVegoAAAAAKqmq0gMAAAAAAAAAAABsC2IpAAAAAAAAAACgEMRSAAAAAAAAAABAIdRUegAAAAAAAAAAqITLb5tXtrXOmja5bGsBsPW4sxQAAAAAAAAAAFAIYikAAAAAAAAAAKAQPIYPWM+9Ty6r9AhsgfnL2rOwrzU9XUtT2shtYt3uFQAAAAAAAAD+yp2lAAAAAAAAAACAQnBnKQAAAAAAAACGjLlPPJdla5and/TKsq1ZN7ohB7/7hLKtB8DwJZYCAAAAAAAAYMhYvaY3tb2dqe99PrUr5g16ve7qMenOLmWYDICRQCwFAAAAAAAAwJDSUL06E2pWZZfSskGt09HVk5a6pEcsBcD/EEsBAAAAAAAAMOTU11bnxCMPHNQa19w6Oy19ZRoIgBGhqtIDAAAAAAAAAAAAbAtiKQAAAAAAAAAAoBDEUgAAAAAAAAAAQCGIpQAAAAAAAAAAgEKoqfQAAAAAAAAAAMPFvU8u26zzlrf3pCq9aarr2+z3AABbnztLAQAAAAAAAAAAheDOUgAAAAAAAADAsDFjxoy0tbVt8Pqsx5du8VpPPHJ/tu99Pg1N3eUYDRgGxFIAAAAAAAAAwLDR1taW1sULkq6W9V6vXdG6xWtVdy5NSp3pXVuu6YChTiwFAAAAAAAAAAwvXS0ptTyTMQ11/S/tUmrf4mUWrGlNVU2SjCrfbJTF5bfNq/QIjFBiKQAAAAAAAABg2BnTUJezjz+0/9f3Prlsi9d4+PFnkvSWcSpgqBNLAYxgamsAAAAAAAAA+CuxFAAAAAAAAADACDf3ieeybM3y9I5eWdZ160Y35OB3n1DWNWFrEksBAAAAAAAAAIxwq9f0pra3M/W9z6d2RXmeUNNdPSbd2aUsa8G2IpYCAAAAAAAAACiAhurVmVCzKruUlg16rY6unrTUJT1iKYYZsRQAAAAAAAAAQEHU11bnxCMPHPQ619w6Oy19ZRgItjGxFMAwtrWeK0zlebYzAAAAAAAAQPmJpQCGsa3xXGEq68knn0pXb3VK9WPT3dlR6XEoMxEcAAAAAGx79z65+Y+aWt7ek6r0pqmub4veBwAMH2IpgGGunM8VpvIWdC1PQ2lU6nt7BHAjTHf1mHR7ZjcAAAAAAABUlFgKYAQo13OFqbyHH38mVRHAjTQdXT1pqUt6xFIAAAAAAABQUWIpABiCBHAjyzW3zk5LX6WnAAAAAAAAAKoqPQAAAAAAAAAAAMC2IJYCAAAAAAAAAAAKwWP4AAAAAKAgLr9tXqVHKKtZjy9N7YrW7FJqz71PLqv0OAAAAMAw4M5SAAAAAAAAAABAIYilAAAAAAAAAACAQvAYPgAAAAAAAAAYpHI+9vqsaZPLthYA6xNLAQAAAAAAADAizX3iuSxbszy9o1eWdd260Q05+N0nlHVNALYNsRQAAAAAAAAAI9LqNb2p7e1Mfe/zqV1Rnjs/dVePSXd2KctaAGx7YikAAAAAAAAARqyG6tWZULMqu5SWDXqtjq6etNQlPWIpgGFLLAUAAAAAAADAiFZfW50Tjzxw0Otcc+vstPSVYSAAKqaq0gMAAAAAAAAAAABsC2IpAAAAAAAAAACgEMRSAAAAAAAAAABAIdRUegAAAAAAAAAAoHguv23egN436/GlqV3Rml1K7bn3yWVlngoY6dxZCgAAAAAAAAAAKASxFAAAAAAAAAAAUAhiKQAAAAAAAAAAoBBqKj0AAAAAAAAADCX3PrnsZc9Z3t6TqvSmqa7vZc8/eM8dyjUaAACD5M5SAAAAAAAAAABAIYilAAAAAAAAAACAQvAYPgAAAKCsNueRJUPdljxSBQAAAAAYPtxZCgAAAAAAAAAAKAR3lgIAAAAAAACAbezeG3+a7s6OjR7rmzNhwOs2NTWlubl5wO8H1vdSf1cHo250Qw5+9wllX5eXJ5YCAAAAAAAAgG2su7Mjfa0LU9e7aoNjrU8vGtii9eOSTBrcYMB6Xurv6oDXrB6T7uxStvXYMmIpAAAAAAAAAKiAut5VGde9MA31teu9PnZ11xavtaqjO33lGgxYz6b+rg5ER1dPWuqSHrFUxYilAAAAAAAAAKBCGuprc+KRB6732sF77rDF61w28+60lmsoYAMb+7s6ENfcOjstysaKqqr0AAAAAAAAAAAAANuCO0sBAAAAUAhzn3guy9YsT+/olRsc65szYdsPtJU1NTWlubm50mMAAMCI8lL7ii31xCP3Z/ve59PQ1D34wQDYbGIp+P/bu/cgO+v6fuCfs/fdZLMEGiAhEJuxP1GwQx2oFVJNxeJktFdGK8wwamq1ojMinY6jMxWL44z1j2indHDwApnihaKVaXFG7FTMxEnaYAuk8ZKEBJIgAXLb3exms5ez398fmG0iuW2+z7nseV6vmZ0Je5798Hk4X875vrOf8zwAAABAKUxOVaOrOhY91f3RdWjbCY8N73q+QV3VSM9ARFzW6C4AAKDlnC5XzFb72L6IylhUpwtqDoCzYlgKAAAAgNLoa5+MRR2HY3HlwAnfXzB5tEEdFe/wkfFIjW4CAABa2KlyxWztnhqOto6IiM5C+gLg7BiWAgAAAKBUerra46Ybrj7he29YfkGDuinemgd/FMONbgIAAFrcyXLFbD2+dU9EVItpCICz1tboBgAAAAAAAAAAAOrBlaUAAAAAAAAAYI7b+NPdMVzdFz3PHCys5vz58+P9739/YfUAmoFhKQAAAAAAAACY4yYmp6JSHY2uI8/H8K7/zS/YMxARl+XXAWgyhqUAAABouI07DzS6hSwHRyeiLaoxvzud9lzesPyCOnYFAAAAlE1PZTIGpgdjweSLWXUOHxmPVFBPAM3GsBQAAAAAAAAAtIi+7o64/R0rsmqsefBHMVxQPwDNpq3RDQAAAAAAAAAAANSDK0sBAAAAUHpz/Xagx3vmwGjsTcMxcXRfVP59W6PbAQAAAGgqriwFAAAAAAAAAACUgitLAQAAAAAAMOed6UqRB0cnoi2qMb87tdRVJQEAmB3DUgAAAAAANI2N330gxseONLqNuuru7YuP/v7fNLoNAACAUjAsBQBQY1t2PBcHpg5GtXew0a1QA36pAQAAxRofOxJpeG90Vw83upW6GG/vj/FY3Og2AAAASsOwFABAjU1OVaOrOhY91f3RdWhbo9uhQH6pAQBAvZTpQxg7Nm+KhdX9sbhrNC6+ZGmj26mpI0cnYqg7YkKuAACAQhV5xdpjGaVv/ngh9Wg8w1IAAHXQ1z4ZizoOx+LKgUa3QkH8UgMAgHoq04cw2sf2RVTGoquvEjfdcHWj26mpb3z/xzGUGt0FAAC0niKvWHsso1SnC2iMpmBYCgCgTnq62lv+L/rLxC81AACot7J8CGP31HC0dUREdDa6FQAAYA7rrh6OgfG90dfTlVVHRmk9hqUAAAAAoIWU8XZtZboVQhk+hPH41j0RUW10GwAAQAvo6+nKzlAySusxLAUAAAAALaSMt2tzKwQAAADgbBmWAgAAAIAW43ZtANBcNu4883vywdGJaItqzO9Opz3+DcsvKLI1AIDSMSwFAAAAAC3I7doAAAAAXq6t0Q0AAAAAAAAAAADUg2EpAAAAAAAAAACgFNyGDwAAAAAAAABa2MadB2Z1/DMHRmNvGo6Jo/ui8u/bTnjso7///4psDaDuDEsBAAAAAAAAADVzzz33xMjIyMu+v27rvnOqt2PzplhY3R9988dzWwNKyLAUAAAAAAAAAFAzIyMjMfzC7oijQyd8v+vQ8DnVax/bF1EZi+p0Ed0BZWNYCgAAAAAAAACoraNDURnaE/193TPfWlwZPadSu6eGo60jIqKzmN6AUjEsBQAAAAAAAABNZOPOA7P+mYOjE9EW1Zjfnc7p5+uhv687bn/Hipl/Ptc+H9+6JyKqBXUFlE1boxsAAAAAAAAAAACoB1eWAgAAAAAAAACAOtiy47k4MHUwqr2Dhdbt7u2Lj/7+3xRas1UZlgIAAAAAAAAAgDqYnKpGV3Useqr7o+vQtkJqjrf3x3gsLqRWGRiWAgAAAAAAAACAOulrn4xFHYdjceVAdq0jRydiqDtiwrDUWTMsBQAAAAAAAAAAddTT1R433XB1dp1vfP/HMZQKaKhE2hrdAAAAAAAAAAAAQD24shQAAAAAAAC0kI0/3R3D1X3R88zBQuvOnz8/3v/+9xdaEwAiIjZ+94EYHztSWL0dmzfFwur+6Js/XlhNWodhKbJt3Jl/D01q7+DoRLRFNeZ3J88ZAAAAAAC0sInJqahUR6PryPMxvOt/iynaMxARlxVTCwB+xfjYkUjDe6O7eriQeu1j+yIqY1GdLqQcLcawFAAAAAAAALSYnspkDEwPxoLJF7NrHT4yHqmAngDgdLqrh2NgfG/09XRl19o9NRxtHRERndm1aD2GpQAAAAAAAKAF9XV3xO3vWJFdZ82DP4rhAvoBgDPp6+mKm264OrvO41v3REQ1vyFaUlujGwAAAAAAAAAAAKgHw1IAAAAAAAAAAEApuA0fAAAAAAAAzBEbdx444zEHRyeiLaoxvzud9vg3LL+gyNYAAOYEV5YCAAAAAAAAAABKwbAUAAAAAAAAAABQCoalAAAAAAAAAACAUuhodAMAAAAAAAAAwNzw+X/fNuufWbd1X3QdGo7FldHYuPNADboCOHuGpQAAAAAAAAAAmJUtO56LA1MHo9o7mF1rx+ZNsbC6P/rmj+c3BmdgWAoAAAAAAAAAgFmZnKpGV3Useqr7o+vQ7K84drz2sX0RlbGoThfUHJyGYSkAAAAAAAAAAGatr30yFnUcjsWVvNsr7p4ajraOiIjOQvqC0zEsBQAAAAAAAADAOenpao+bbrg6q8bjW/dERLWYhuAM2hrdAAAAAAAAAAAAQD24shQAAAAAAABQF/fcc0+MjIwUVm/dunXRdeT5WNozVlhNAKC1GZYCAAAAAAAA6mJkZCSGX9gdcXSokHpHD+6NShqNamch5QCAEjAsBQAAAAAAANTP0aGoDO2J/r7uAmoN/vI3nqalAICzY1gKAAAAAAAAqKv+vu64/R0rTvn4xp0HzqrOdzb8PI5OVmN+dzrrnwEAyq2t0Q0AAAAAAAAAAADUg2EpAAAAAAAAAACgFAxLAQAAAAAAAAAApWBYCgAAAAAAAAAAKAXDUgAAAAAAAAAAQCkYlgIAAAAAAAAAAErBsBQAAAAAAAAAAFAKhqUAAAAAAAAAAIBSMCwFAAAAAAAAAACUQkejGwAAAICy2LjzQKNbyHJwdCLaohrzu9OcPxcAAAAAoJxcWQoAAAAAAAAAACgFV5YCAAAAAAAATmnjT3fHcHVf9DxzMLvWunXrouvI87G0Z6yAzgAAZs+wFAAAAAAAAHBKE5NTUamORteR52N41/9m1Tp6cG9U0mhUOwtqDgBglgxLAQAAAAAAAKfVU5mMgenBWDD5Yl6ho4O//A2laSkAoDEMSwEAAAAAAABn1NfdEbe/Y0VWjXVP7oyYrhbUEQDA7LU1ugEAAAAAAAAAAIB6cGUpAAAAAAAAKKGNOw+c1XEHRyeiLaoxvzud9mfesPyColoDGmzLjufiwNTBqPYOFlJvx+ZNsbC6P/rmjxdSDyCHYSkAAAAAAAAAYMbkVDW6qmPRU90fXYe2ZddrH9sXURmL6nQBzQFkMiwFAAAAAAAAAJygr30yFnUcjsWVs7sK3ensnhqOto6IiM7sWgC5DEsBAAAAAAAAAC/T09UeN91wdXadx7fuiYhqfkMABWhrdAMAAAAAAAAAAAD1YFgKAAAAAAAAAAAoBcNSAAAAAAAAAABAKRiWAgAAAAAAAAAASqGj0Q0AAAAAAAAAc9/GnQfOeMzB0Yloi2rM705ndTwAQNFcWQoAAAAAAAAAACgFw1IAAAAAAAAAAEApGJYCAAAAAAAAAABKwbAUAAAAAAAAAABQCoalAAAAAAAAAACAUjAsBQAAAAAAAAAAlEJHoxv4VVNTUzN/3rt3bwM7KbfBwcEYOTwa01Nj8ey+odMe++LBw3XqihzjE1PRFtU4Ml7xnLUQz2vr8Zy2Js9raxo5Mh5H0pGY7ByOZ599ttHtzDh+D3383rps5IrmUKZcUZbX+rKcZ4RzbVVlOdeynGeEc21VZTpXuaK5yRXNochcUfTrSzPX01tz1NNbc9TTW+NrNXs9vTVHPb2dO7li9ioppdToJo732GOPxW//9m83ug0AAJjzNm3aFNdcc02j22gIuQIAAIohV8gVAACQq9lyRdPdhu+FF15odAsAANASyry3LvO5AwBAkcq8ty7zuQMAQJGabW/ddLfhu/zyy2f+vGHDhrj00ksb2A1z0d69e2c+7bNp06ZYvHhxgztirrGGyGH9kMsaIteePXvi2muvjYgT99ZlI1eQy+sxuawhclg/5LKGyCVXvESuIJfXY3JZQ+SwfshlDZGrmXNF0w1L9fT0zPz50ksvjaVLlzawG+a6xYsXW0NksYbIYf2Qyxoi1/F767KRKyiS12NyWUPksH7IZQ2RS654iVxBLq/H5LKGyGH9kMsaIlez5Yqmuw0fAAAAAAAAAABALRiWAgAAAAAAAAAASsGwFAAAAAAAAAAAUAqGpQAAAAAAAAAAgFIwLAUAAAAAAAAAAJSCYSkAAAAAAAAAAKAUDEsBAAAAAAAAAAClUEkppUY3AQAAAAAAAAAAUGuuLAUAAAAAAAAAAJSCYSkAAAAAAAAAAKAUDEsBAAAAAAAAAAClYFgKAAAAAAAAAAAoBcNSAAAAAAAAAABAKRiWAgAAAAAAAAAASsGwFAAAAAAAAAAAUAqGpQAAAAAAAAAAgFIwLAUAAAAAAAAAAJRCXYalNm3aFO9617vikksuie7u7rjooovi7W9/e3zve98755rf/OY34/rrr49FixZFT09PLFu2LN773vfGz372swI7pxnUYv18+9vfjje+8Y2xYMGC6O7ujmXLlsW73/3u+OlPf1pg5zSLWqyh4z3yyCNRqVRmvn74wx8WUpfmUMv1s23btnj1q18dlUolrrnmmvjFL35RQMc0m6LX0Pbt2+NDH/pQvOY1r4l58+ZFV1dXLFmyJP7gD/4gHnrooWKbp2l87GMfO+G95otf/GJWvbm6l5YryCFXkEuuIIdcQS65giLIFS+RK8ghV5BLriCHXEEuuYIitEyuSDV23333pfb29hQRKSJSb2/vzJ8jIv3t3/7trGu+733vm/n5SqWSenp6Zv65p6cnPfLIIzU4Exqh6PUzPj6eVq1adUKNzs7OmT/39fWldevW1ehsaIRavAb9qmuuuSa1tbWlSqWSIiI9+uij+Y3TFGq5fn7yk5+kCy+8MEVE+tM//dM0NjZWYOc0i6LX0He/+93U3d19wj6oq6vrhJrvfe97a3Q2NML09HS69dZbU0SkhQsXzqyhu++++5xrztW9tFxBDrmCXHIFOeQKcskV5JIr/o9cQQ65glxyBTnkCnLJFeRqtVxR02Gp7du3z/wP8fGPfzy9+OKLKaWUBgcH05o1a1JHR0eqVCqz2uytXbs2RUTq7+9P9913Xzpy5EhKKaVdu3al1atXp4hI559/fhocHKzJOVE/tVg/73nPe1JEpEsuuSR961vfSqOjo2lycjJt3LgxvfrVr04RkZYvX56q1WqtTos6qsUa+lXf+c53UkSk2267LS1btkz4aCG1XD8vvPBCWrp0aYqI9O53v9trTosqeg2NjY3NBNbXvOY16Qc/+EE6evRomp6eTk8//XT68z//85kN5EMPPVTLU6NOpqamZvYuF154YXriiSdm3mvONXzM1b20XEEOuYJccgU55ApyyRXkkiv+j1xBDrmCXHIFOeQKcskV5GrFXFHTYakPfvCDKSLS+973vpM+/ulPfzpFRFq1atVZ1zy2Qbz//vtP+vhb3vKWFBHp7/7u786pZ5pH0etn69atqa2tLXV2dqYnnnjiZY//+Mc/nnnRXr9+fVbvNIdavAYdb3p6Or32ta9Nvb29ae/evcJHi6nl+rn++utTRKQ//MM/FDxaWNFr6Ac/+MHM+9STTz75ssenp6fT6173uplQy9w2MTGR3vnOd878penPf/7zlFLKDh9zdS8tV5BDriCXXEEOuYJccgU55IoTyRXkkCvIJVeQQ64gl1xBjlbNFTUdllqyZEmKiLRly5aTPn7o0KHU0dGROjs708jIyBnrbd26dWZS7VQv1g8//HCKiHTttddm9U7jFb1+du/enT7ykY+k1atXn/KYiy++OEVEuuuuu865b5pH0WvoV33961+f+ZRGSkn4aDG1Wj9f/epXU0Skyy67LA0PDxfVLk2o6DV07DWnv7//lMd8+MMfThGR3vrWt55z3zSHzZs3p97e3rR8+fL09NNPz3w/J3zM5b20XEEOuYJccgU55ApyyRXkkCtOJFeQQ64gl1xBDrmCXHIFOVo1V7RFjezYsSOee+65WLhwYVxxxRUnPea8886LK6+8MiYnJ+M///M/z1hz/fr1ERFx3XXXRVvbyVtfsWJFRERs2rQpxsfHz7F7Gq0W6+fSSy+NL3zhC/GVr3zllMcsXbo0IiKGhobOrXGaRi3W0PGq1Wrccccd0dvbGx/72MeKaJkmUqv1Mzk5GXfccUdERHzuc5+L/v7+wnqmudRiDV188cVn/e+fzbE0p9e+9rXx0EMPxfr16+MVr3hFITXn6l5ariCHXEEuuYIccgW55ApyyRX/R64gh1xBLrmCHHIFueQKcrVqrqjZsNT27dsj4qUN3+kce/ypp54qpObAwED09/fH1NRUPPPMM2fZLc2mFuvnbExMTETES28IzG21XkP33XdfbN++PT7wgQ94k29BtVo/X//612PPnj1x3XXXxZ/92Z/lNUlTq8Uaet3rXhcLFiyIw4cPx+bNm1/2eEopNmzYEBERK1eunGXHNKMbbrghlixZUli9ubqXlivIIVeQS64gh1xBLrmCIsgVL5EryCFXkEuuIIdcQS65giK0Yq6o2bDU7t27IyKit7f3tMf19fWdcHy9a9KcGvFcT0xMxM6dOyMi4td//dez69FYtVxDExMTceedd/qURgur1fpZu3ZtRER89KMfjQ0bNsRb3/rWuOiii+K8886La6+9Nr761a9GSimjc5pFLdbQwMBAfOYzn4mIiJtuuinWrVsXExMTkVKKXbt2xV/8xV/E//zP/8Tv/M7vxC233JJ5BrSiubqXlivIIVeQS64gh1xBLrmCZjRX99JyBTnkCnLJFeSQK8glV9CMmmEv3VF4xV8aGRkp/Pha1KQ5NeK5/td//dcYGRmJBQsWxBvf+MbsejRWLdfQPffcE7t3747bbrvNpzRaVC3Wz4EDB2LdunUxb9686OjoiJUrV8bk5OTM4xs3boyNGzfGo48+Gv/0T/80655pLrV6Dfrwhz8cF110UXz2s5+NlStXRqVSic7OzpiYmIjzzz8//vqv/zo+9alPRXt7+7m0TYubq3tpuYIccgW55ApyyBXkkitoRnN1Ly1XkEOuIJdcQQ65glxyBc2oGfbSNbuy1NjY2KyOP3LkSENq0pzq/VwfPXo0Pv7xj0dExEc+8pGYN29eVj0ar1ZraGxsLD7zmc/4lEaLq8X6+a//+q+Ynp6OK6+8Mj75yU/GF7/4xdi3b1+Mj4/Hk08+GTfeeGNERNx///3xta997Zz6pnnU8n1samoqBgYGIuKlS9keuyR7b29v9Pb2zvrfTXnM1b20XEEOuYJccgU55ApyyRU0o7m6l5YryCFXkEuuIIdcQS65gmbUDHvpmg1LzfayfJVKpSE1aU71fq7/8i//Mp566qm48sor4xOf+ERWLZpDrdbQP/zDP8Tzzz/v3t8trhbr54knnoiIl0LIZz/72Vi9enX82q/9WnR1dcVv/uZvxoMPPjjzKbG777571j3TXGq1D1q9enXcfPPN8eyzz8a3vvWt+MUvfhHDw8Px2GOPxXXXXRd33nlnvP71r5+5TDscb67upeUKcsgV5JIryCFXkEuuoBnN1b20XEEOuYJccgU55ApyyRU0o2bYS9dsWOrYvQPP1pnuRVirmjSnej7Xd9xxR6xduzYWLlwY//Iv/xI9PT3nXIvmUYs1dPjw4fjc5z7nUxolUIv1s3///oiIeMUrXhGrVq162eOVSiVuvfXWiIh47LHHTrjkLXNPLdbQN7/5zbj33ntj4cKFsX79+rjxxhtjyZIl0d/fH1dffXU88MAD8c53vjN27Ngxs5bgeHN1Ly1XkEOuIJdcQQ65glxyBc1oru6l5QpyyBXkkivIIVeQS66gGTXDXrpmw1Lz588v/Pha1KQ51eu5XrNmTdx5553R29sb//Zv/xa/8Ru/cU51aD61WENr1qyJAwcO+JRGCdRi/QwPD0dExJVXXnnKY6644oqIiJiYmIgXX3xxVj3QXGqxhh544IGIiLjpppvioosuOukxt99+e0REfP/7349Dhw7Nqgda31zdS8sV5JAryCVXkEOuIJdcQTOaq3tpuYIccgW55ApyyBXkkitoRs2wl67ZsNRll10WEWe+1+CxewseO77eNWlO9Xiu77nnnvirv/qr6OjoiAcffDCuu+662TdK0yp6DR08eDDWrFnjUxolUYvXoGMTz/PmzTvlMce/0R89evSMNWletVhD27Zti4g47V+UvepVr4qIly5f+tRTT51Vr5THXN1LyxXkkCvIJVeQQ64gl1xBM5qre2m5ghxyBbnkCnLIFeSSK2hGzbCX7ii84i+98pWvjIiIPXv2nPa4Y48fOz635tDQUBw+fDg6Ojpi2bJlZ9suTaYW6+d4X/va1+KDH/xgVCqVWLt2bbztbW87t0ZpWkWvoc2bN8fw8HC0t7fH5ZdfftJjjk3iv/3tb4+Ojo5YsWJFPPzww7NtnSZQi9egY5P1p3vTP/aGHxFx3nnnnbEmzasWa6harUZERFvbqWfdj79n87Hj4Zi5upeWK8ghV5BLriCHXEEuuYJmNFf30nIFOeQKcskV5JAryCVX0IyaYS9dsytLvfKVr4zFixfHoUOH4ic/+clJjxkcHIwtW7ZEZ2dnvP71rz9jzd/93d+NiIgNGzbE9PT0SY/50Y9+FBER11xzjXs5z2G1WD/HPPTQQ/Ge97wnpqen46677oqbb765qLZpIrVaQ9VqNYaGhk76lVKKiIjR0dEYGhqKkZGRws6H+qrF+rnqqqsiIuJnP/vZKY859tiiRYviggsumH3jNI1arKHFixdHRMT27dtPecyxT3NERCxZsmSWXdPq5upeWq4gh1xBLrmCHHIFueQKmtFc3UvLFeSQK8glV5BDriCXXEEzaoa9dM2GpSIi/uiP/igiIr7whS+c9PG77rorpqam4i1vectZ3WPwVa96VVx++eXxwgsvxDe+8Y2THvP3f//3ERHxx3/8x+fUM82j6PUT8dI9Ud/1rnfF1NRUfPrTn45bb721qHZpQkWuoZUrV0ZK6bRfxyZaH3300UgpxQ9/+MMiT4c6K/o16Prrr4/+/v7Yvn17rFu37qTHfOlLX4qIiFWrVp1b0zSVotfQ7/3e70VExD//8z/HwYMHT3rMP/7jP0ZExPLly13en5eZy3tpuYIccgW55ApyyBXkkitoNnN5Ly1XkEOuIJdcQQ65glxyBc2mKfbSqYa2bduWurq6UkSkT3ziE2nfvn0ppZSGhobS5z//+dTR0ZEqlUpat27dzM8MDg6m3/qt30o9PT3plltueVnNtWvXpohI/f39ae3atWlsbCyllNLu3bvT6tWrU0Sk888/Pw0ODtby1KiDotfP448/nvr6+lJEpNtvv72u50Jj1OI16HSWLVuWIiI9+uijRZ4GDVKL9fPJT34yRUS68MIL07e//e00Pj6eUjrxPayzszNt3ry5PidJTRW9hp5//vk0MDCQIiJdddVV6T/+4z/S6Ohoqlar6amnnkof+MAHUkSkiEhf+tKX6nqu1M+x95q77777pI+36l5ariCHXEEuuYIccgW55ApqQa6QK5g9uYJccgU55ApyyRXUwlzPFTUdlkoppXvvvTe1t7fP/M9wbPN37OtTn/rUCcevX79+5rGOjo6T1jz2HyYiUqVSSb29vTP/3NPTk773ve/V+rSokyLXz7333jvz2MDAwBm/du3aVc9TpUZq8Rp0KsJH6yl6/YyPj6e3ve1tM8e0t7efULOtrS195StfqdfpUQdFr6FHHnkk9ff3n7AP6uzsPKHmhz70oXqdHjX25je/ObW3t5/wdfzrxfHff/Ob35xSau29tFxBDrmCXHIFOeQKcskV5JArTiRXkEOuIJdcQQ65glxyBTlaMVfUfFgqpZQ2bNiQbrzxxrR48eLU2dmZFi1alFatWpUefvjhlx07ODiYrrrqqtTd3X3aKen7778/rVy5Mp1//vmpq6srXXrppemWW25JW7ZsqeWp0ABFrZ/jw8fZfD399NN1OkNqrRavQScjfLSmotdPtVpNX/7yl9OKFSvSwMBA6urqSkuXLk0333xz+u///u9anw4NUPQaeuaZZ9Jtt92WrrjiijRv3rzU2dmZlixZkv7kT/7EX8C2mDe96U1nvW9505velFJq/b20XEEOuYJccgU55ApyyRWcK7ni5eQKcsgV5JIryCFXkEuu4Fy1Yq6opJRSAAAAAAAAAAAAtLi2RjcAAAAAAAAAAABQD4alAAAAAAAAAACAUjAsBQAAAAAAAAAAlIJhKQAAAAAAAAAAoBQMSwEAAAAAAAAAAKVgWAoAAAAAAAAAACgFw1IAAAAAAAAAAEApGJYCAAAAAAAAAABKwbAUAAAAAAAAAABQCoalAAAAAAAAAACAUjAsBQAAAAAAAAAAlIJhKQAAAAAAAAAAoBQMSwEAAAAAAAAAAKVgWAoAAAAAAAAAACgFw1IAAAAAAAAAAEApGJYCAAAAAAAAAABKwbAUAAAAAAAAAABQCoalAAAAAAAAAACAUjAsBQAAAAAAAAAAlIJhKQAAAAAAAAAAoBQMSwEAAAAAAAAAAKVgWAoAAAAAAAAAACgFw1IAAAAAAAAAAEApGJYCAAAAAAAAAABKwbAUAAAAAAAAAABQCv8fPf8IUwQX4uYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1200x400 with 3 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 390,
       "width": 1189
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(1, 3, figsize=(12, 4))\n",
    "ax[0].hist(pscores[z == 1], bins=5, alpha=0.5, density=True)\n",
    "ax[0].hist(pscores[z == 0], bins=5, alpha=0.5, density=True, edgecolor=\"black\")\n",
    "ax[1].hist(pscores[z == 1], bins=10, alpha=0.5, density=True)\n",
    "ax[1].hist(pscores[z == 0], bins=10, alpha=0.5, density=True, edgecolor=\"black\")\n",
    "ax[2].hist(pscores[z == 1], bins=30, alpha=0.5, density=True, label=\"treated\")\n",
    "ax[2].hist(\n",
    "    pscores[z == 0],\n",
    "    bins=30,\n",
    "    alpha=0.5,\n",
    "    density=True,\n",
    "    label=\"control\",\n",
    "    edgecolor=\"black\",\n",
    ")\n",
    "ax[2].legend()\n",
    "\n",
    "# Hide the y-axis\n",
    "ax[0].set_ylim(0, 4.5), ax[1].set_ylim(0, 4.5), ax[2].set_ylim(0, 4.5)\n",
    "ax[0].set_xlim(0, 1), ax[1].set_xlim(0, 1), ax[2].set_xlim(0, 1)\n",
    "ax[0].set_yticks([]), ax[1].set_yticks([]), ax[2].set_yticks([])\n",
    "ax[0].set_yticklabels([]), ax[1].set_yticklabels([]), ax[2].set_yticklabels([])\n",
    "ax[0].set_title(\"5 bins\"), ax[1].set_title(\"10 bins\"), ax[2].set_title(\"30 bins\")\n",
    "f.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alal/Desktop/code/00_causal/CausalinferencePy/causalinference/core/summary.py:102: RuntimeWarning: invalid value encountered in divide\n",
      "  return (mean_t - mean_c) / np.sqrt((sd_c**2 + sd_t**2) / 2)\n",
      "/home/alal/Desktop/code/00_causal/CausalinferencePy/causalinference/estimators/ols.py:21: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  olscoef = np.linalg.lstsq(Z, Y)[0]\n",
      "/home/alal/Desktop/code/00_causal/CausalinferencePy/causalinference/core/summary.py:28: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  self._dict[\"Y_t_sd\"] = np.sqrt(data[\"Y_t\"].var(ddof=1))\n",
      "/home/alal/anaconda3/envs/metrics/lib/python3.11/site-packages/numpy/core/_methods.py:261: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/alal/Desktop/code/00_causal/CausalinferencePy/causalinference/core/summary.py:33: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  self._dict[\"X_t_sd\"] = np.sqrt(data[\"X_t\"].var(0, ddof=1))\n",
      "/home/alal/anaconda3/envs/metrics/lib/python3.11/site-packages/numpy/core/_methods.py:258: RuntimeWarning: invalid value encountered in divide\n",
      "  ret = um.true_divide(\n",
      "/home/alal/Desktop/code/00_causal/CausalinferencePy/causalinference/core/summary.py:27: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  self._dict[\"Y_c_sd\"] = np.sqrt(data[\"Y_c\"].var(ddof=1))\n",
      "/home/alal/Desktop/code/00_causal/CausalinferencePy/causalinference/core/summary.py:32: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  self._dict[\"X_c_sd\"] = np.sqrt(data[\"X_c\"].var(0, ddof=1))\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>5</th>\n",
       "      <th>10</th>\n",
       "      <th>20</th>\n",
       "      <th>50</th>\n",
       "      <th>80</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>est</th>\n",
       "      <td>-0.116096</td>\n",
       "      <td>-0.177604</td>\n",
       "      <td>-0.199688</td>\n",
       "      <td>-0.264742</td>\n",
       "      <td>-0.203770</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>se</th>\n",
       "      <td>0.281897</td>\n",
       "      <td>0.279457</td>\n",
       "      <td>0.272488</td>\n",
       "      <td>0.256721</td>\n",
       "      <td>0.244644</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           5         10        20        50        80\n",
       "est -0.116096 -0.177604 -0.199688 -0.264742 -0.203770\n",
       "se   0.281897  0.279457  0.272488  0.256721  0.244644"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = CausalModel(y, z, X)\n",
    "m.est_propensity()\n",
    "m.cutoff = 0\n",
    "res = []\n",
    "for nn in [5, 10, 20, 50, 80]:\n",
    "    m.blocks = nn\n",
    "    m.stratify()\n",
    "    m.est_via_blocking(adj=0)\n",
    "    ests = m.estimates[\"blocking\"]\n",
    "    res.append(np.array([ests[\"ate\"], ests[\"ate_se\"]]))\n",
    "\n",
    "pd.DataFrame(np.c_[res].T, index=[\"est\", \"se\"], columns=[5, 10, 20, 50, 80])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## weighting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ate': -1.5162837769469517, 'ate_se': nan}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.reset()\n",
    "m.est_propensity()\n",
    "m.cutoff = 0\n",
    "m.trim()\n",
    "m.est_via_weighting(estimand=\"ate\")\n",
    "m.estimates[\"weighting\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ate': -0.15566888253380995, 'ate_se': nan}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.reset()\n",
    "m.est_propensity()\n",
    "m.cutoff = 0\n",
    "m.trim()\n",
    "m.est_via_weighting(estimand=\"ate\", hajekize=True)\n",
    "m.estimates[\"weighting\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ate': -1.5162837769469517, 'ate_se': nan}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.reset()\n",
    "m.est_propensity()\n",
    "m.cutoff = 0.01\n",
    "m.trim()\n",
    "m.est_via_weighting(estimand=\"ate\")\n",
    "m.estimates[\"weighting\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ate': -0.15566888253380995, 'ate_se': nan}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.reset()\n",
    "m.est_propensity()\n",
    "m.cutoff = 0.01\n",
    "m.trim()\n",
    "m.est_via_weighting(estimand=\"ate\", hajekize=True)\n",
    "m.estimates[\"weighting\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ate': -1.3034023808081816, 'ate_se': nan}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.reset()\n",
    "m.est_propensity()\n",
    "m.cutoff = 0.05\n",
    "m.trim()\n",
    "m.est_via_weighting(estimand=\"ate\")\n",
    "m.estimates[\"weighting\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ate': -0.10280308535143234, 'ate_se': nan}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.reset()\n",
    "m.est_propensity()\n",
    "m.cutoff = 0.05\n",
    "m.trim()\n",
    "m.est_via_weighting(estimand=\"ate\", hajekize=True)\n",
    "m.estimates[\"weighting\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ate': -0.0859162644736197, 'ate_se': nan}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.reset()\n",
    "m.est_propensity()\n",
    "m.cutoff = 0.1\n",
    "m.trim()\n",
    "m.est_via_weighting(estimand=\"ate\")\n",
    "m.estimates[\"weighting\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ate': 0.15180636527279034, 'ate_se': nan}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.reset()\n",
    "m.est_propensity()\n",
    "m.cutoff = 0.1\n",
    "m.trim()\n",
    "m.est_via_weighting(estimand=\"ate\", hajekize=True)\n",
    "m.estimates[\"weighting\"]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "metrics",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}