{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python: Sample Selection Models\n",
    "\n",
    "In this example, we illustrate how the [DoubleML](https://docs.doubleml.org/stable/index.html) package can be used to estimate the average treatment effect (ATE) under sample selection or outcome attrition. The estimation is based on a simulated DGP from Appendix E of [Bia, Huber and Lafférs (2023)](https://doi.org/10.1080/07350015.2023.2271071).  \n",
    "\n",
    "Consider the following DGP:\n",
    "$$\n",
    "\\begin{align*}\n",
    "Y_i &= \\theta_0 D_i + X_i'\\beta_0 + \\varepsilon_i,\\\\\n",
    "S_i &= \\mathbb{1}\\{D_i + \\gamma_0 Z_i + X_i'\\beta_0 + \\upsilon_i > 0\\}, \\\\\n",
    "D_i &= \\mathbb{1}\\{X_i'\\beta_0 + \\xi_i > 0\\}\n",
    "\\end{align*}\n",
    "$$\n",
    "where $Y_i$ is observed if $S_i=1$\n",
    "with\n",
    "$$X_i \\sim N(0, \\sigma^2_X), \\quad Z_i \\sim N(0, 1), \\quad (\\varepsilon,_i \\nu_i) \\sim N(0, \\sigma^2_{\\varepsilon, \\nu}), \\quad \\xi_i \\sim N(0, 1).$$\n",
    "\n",
    "Let $D_i\\in\\{0,1\\}$ denote the treatment status of unit $i$ and let $Y_{i}$ be the outcome of interest of unit $i$.\n",
    "Using the potential outcome notation, we can write $Y_{i}(d)$ for the potential outcome of unit $i$ and treatment status $d$. Further, let $X_i$ denote a vector of pre-treatment covariates.  \n",
    "\n",
    "## Outcome missing at random (MAR)  \n",
    "Now consider the first setting, in which the outcomes are missing at random (MAR), according to assumptions in [Bia, Huber and Lafférs (2023)](https://doi.org/10.1080/07350015.2023.2271071). \n",
    "Let the covariance matrix $\\sigma^2_X$ be such that $a_{ij} = 0.5^{|i - j|}$, $\\gamma_0 = 0$, $\\sigma^2_{\\varepsilon, \\upsilon} = \\begin{pmatrix} 1 & 0 \\\\  0 & 1 \\end{pmatrix}$ and finally, let the vector of coefficients $\\beta_0$ resemble a quadratic decay of coefficients importance; $\\beta_{0,j} = 0.4/j^2$ for $j = 1, \\ldots, p$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data\n",
    "\n",
    "We will use the implemented data generating process `make_ssm_data` to generate data according to the simulation in Appendix E of [Bia, Huber and Lafférs (2023)](https://doi.org/10.1080/07350015.2023.2271071). The true ATE in this DGP is equal to $\\theta_0=1$ (it can be changed by setting the parameter `theta`). \n",
    "\n",
    "The data generating process `make_ssm_data` by default settings already returns a `DoubleMLData` object (however, it can return a pandas DataFrame or a NumPy array if `return_type` is specified accordingly). In this first setting, we are estimating the ATE under missingness at random, so we set `mar=True`.\n",
    "The selection indicator `S` can be set via `s_col`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================== DoubleMLSSMData Object ==================\n",
      "Selection variable: s\n",
      "Outcome variable: y\n",
      "Treatment variable(s): ['d']\n",
      "Covariates: ['X1', 'X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8', 'X9', 'X10', 'X11', 'X12', 'X13', 'X14', 'X15', 'X16', 'X17', 'X18', 'X19', 'X20', 'X21', 'X22', 'X23', 'X24', 'X25', 'X26', 'X27', 'X28', 'X29', 'X30', 'X31', 'X32', 'X33', 'X34', 'X35', 'X36', 'X37', 'X38', 'X39', 'X40', 'X41', 'X42', 'X43', 'X44', 'X45', 'X46', 'X47', 'X48', 'X49', 'X50', 'X51', 'X52', 'X53', 'X54', 'X55', 'X56', 'X57', 'X58', 'X59', 'X60', 'X61', 'X62', 'X63', 'X64', 'X65', 'X66', 'X67', 'X68', 'X69', 'X70', 'X71', 'X72', 'X73', 'X74', 'X75', 'X76', 'X77', 'X78', 'X79', 'X80', 'X81', 'X82', 'X83', 'X84', 'X85', 'X86', 'X87', 'X88', 'X89', 'X90', 'X91', 'X92', 'X93', 'X94', 'X95', 'X96', 'X97', 'X98', 'X99', 'X100']\n",
      "Instrument variable(s): None\n",
      "No. Observations: 2000\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from doubleml.irm.datasets import make_ssm_data\n",
    "import doubleml as dml\n",
    "\n",
    "np.random.seed(42)\n",
    "n_obs = 2000\n",
    "df = make_ssm_data(n_obs=n_obs, mar=True, return_type='DataFrame')\n",
    "\n",
    "dml_data = dml.DoubleMLSSMData(df, 'y', 'd', s_col='s')\n",
    "print(dml_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimation\n",
    "\n",
    "To estimate the ATE under sample selection, we will use the `DoubleMLSSM` class. \n",
    "\n",
    "As for all `DoubleML` classes, we have to specify learners, which have to be initialized first.\n",
    "Given the simulated quadratic decay of coefficients importance, Lasso regression should be a suitable option (as for propensity scores, this will be a $\\mathcal{l}_1$-penalized Logistic Regression). \n",
    "\n",
    "The learner `ml_g` is used to fit conditional expectations of the outcome $\\mathbb{E}[Y_i|D_i, S_i, X_i]$, whereas the learners `ml_m` and `ml_pi` will be used to estimate the treatment and selection propensity scores $P(D_i=1|X_i)$ and $P(S_i=1|D_i, X_i)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LassoCV, LogisticRegressionCV\n",
    "\n",
    "ml_g = LassoCV()\n",
    "ml_m = LogisticRegressionCV(penalty='l1', solver='liblinear')\n",
    "ml_pi = LogisticRegressionCV(penalty='l1', solver='liblinear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `DoubleMLSSM` class can be used as any other `DoubleML` class. \n",
    "\n",
    "The score is set to `score='missing-at-random'`, since the parameters of the DGP were set to satisfy the assumptions of outcomes missing at random. Further, since the simulation in [Bia, Huber and Lafférs (2023)](https://doi.org/10.1080/07350015.2023.2271071) uses normalization of inverse probability weights, we will apply the same setting by `normalize_ipw=True`.\n",
    "\n",
    "After initialization, we have to call the `fit()` method to estimate the nuisance elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================== DoubleMLSSM Object ==================\n",
      "\n",
      "------------------ Data Summary      ------------------\n",
      "Outcome variable: y\n",
      "Treatment variable(s): ['d']\n",
      "Covariates: ['X1', 'X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8', 'X9', 'X10', 'X11', 'X12', 'X13', 'X14', 'X15', 'X16', 'X17', 'X18', 'X19', 'X20', 'X21', 'X22', 'X23', 'X24', 'X25', 'X26', 'X27', 'X28', 'X29', 'X30', 'X31', 'X32', 'X33', 'X34', 'X35', 'X36', 'X37', 'X38', 'X39', 'X40', 'X41', 'X42', 'X43', 'X44', 'X45', 'X46', 'X47', 'X48', 'X49', 'X50', 'X51', 'X52', 'X53', 'X54', 'X55', 'X56', 'X57', 'X58', 'X59', 'X60', 'X61', 'X62', 'X63', 'X64', 'X65', 'X66', 'X67', 'X68', 'X69', 'X70', 'X71', 'X72', 'X73', 'X74', 'X75', 'X76', 'X77', 'X78', 'X79', 'X80', 'X81', 'X82', 'X83', 'X84', 'X85', 'X86', 'X87', 'X88', 'X89', 'X90', 'X91', 'X92', 'X93', 'X94', 'X95', 'X96', 'X97', 'X98', 'X99', 'X100']\n",
      "Instrument variable(s): None\n",
      "No. Observations: 2000\n",
      "\n",
      "\n",
      "------------------ Score & Algorithm ------------------\n",
      "Score function: missing-at-random\n",
      "\n",
      "------------------ Machine Learner   ------------------\n",
      "Learner ml_g: LassoCV()\n",
      "Learner ml_pi: LogisticRegressionCV(penalty='l1', solver='liblinear')\n",
      "Learner ml_m: LogisticRegressionCV(penalty='l1', solver='liblinear')\n",
      "Out-of-sample Performance:\n",
      "Regression:\n",
      "Learner ml_g_d0 RMSE: [[1.10039862]]\n",
      "Learner ml_g_d1 RMSE: [[1.11071087]]\n",
      "Classification:\n",
      "Learner ml_pi Log Loss: [[0.53791422]]\n",
      "Learner ml_m Log Loss: [[0.63593298]]\n",
      "\n",
      "------------------ Resampling        ------------------\n",
      "No. folds: 5\n",
      "No. repeated sample splits: 1\n",
      "\n",
      "------------------ Fit Summary       ------------------\n",
      "       coef  std err        t          P>|t|    2.5 %    97.5 %\n",
      "d  0.997494  0.03113  32.0428  2.765710e-225  0.93648  1.058508\n"
     ]
    }
   ],
   "source": [
    "dml_ssm = dml.DoubleMLSSM(dml_data, ml_g, ml_m, ml_pi, score='missing-at-random', normalize_ipw=True)\n",
    "dml_ssm.fit()\n",
    "\n",
    "print(dml_ssm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confidence intervals at different levels can be obtained via"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     5.0 %    95.0 %\n",
      "d  0.94629  1.048699\n"
     ]
    }
   ],
   "source": [
    "print(dml_ssm.confint(level=0.90))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ATE estimates distribution\n",
    "\n",
    "Here, we add a small simulation where we generate multiple datasets, estimate the ATE and collect the results (this may take some time). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 0/200\n",
      "Iteration: 20/200\n",
      "Iteration: 40/200\n",
      "Iteration: 60/200\n",
      "Iteration: 80/200\n",
      "Iteration: 100/200\n",
      "Iteration: 120/200\n",
      "Iteration: 140/200\n",
      "Iteration: 160/200\n",
      "Iteration: 180/200\n"
     ]
    }
   ],
   "source": [
    "n_rep = 200\n",
    "ATE = 1.0\n",
    "\n",
    "ATE_estimates = np.full((n_rep), np.nan)\n",
    "\n",
    "np.random.seed(42)\n",
    "for i_rep in range(n_rep):\n",
    "    if (i_rep % int(n_rep/10)) == 0:\n",
    "        print(f'Iteration: {i_rep}/{n_rep}')\n",
    "    dml_data = make_ssm_data(n_obs=n_obs, mar=True)\n",
    "    dml_ssm = dml.DoubleMLSSM(dml_data, ml_g, ml_m, ml_pi, score='missing-at-random', normalize_ipw=True)\n",
    "    dml_ssm.fit()\n",
    "\n",
    "    ATE_estimates[i_rep] = dml_ssm.coef.squeeze()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distribution of the estimates takes the following form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAATeVJREFUeJzt3Xl4VOWhP/DvLJmZLLNkmWSykpCQhIQQwiqiLMoiKG6t9ba2ota2Kv5aS7229HoV2yp6Wy2uqK2C1AU3pHUD2ULYtxB2QvZ939eZzMz5/UFJjRAgYZL3zMz38zzzPGRyZubLSZj5cs573lchSZIEIiIiIhlSig5ARERE1B8WFSIiIpItFhUiIiKSLRYVIiIiki0WFSIiIpItFhUiIiKSLRYVIiIiki0WFSIiIpIttegAV8LpdKKyshJ6vR4KhUJ0HCIiIroMkiShra0NERERUCovfszErYtKZWUloqOjRccgIiKiQSgrK0NUVNRFt3HroqLX6wGc/YsaDAbBaYiIiOhytLa2Ijo6uvdz/GLcuqicO91jMBhYVIiIiNzM5QzbEDqYNjY2FgqF4rzb4sWLRcYiIiIimRB6ROXAgQNwOBy9Xx8/fhxz5szBHXfcITAVERERyYXQomI2m/t8/eyzzyI+Ph4zZswQlIiIiIjkRDZjVGw2G959910sWbKk33NWVqsVVqu19+vW1tbhikdERAJJkgS73d7nKDzJl0qlglqtdsnUIbIpKuvXr0dzczPuueeefrdZvnw5nnrqqeELRUREwtlsNlRVVaGzs1N0FBoAPz8/hIeHQ6PRXNHzKCRJklyU6YrMmzcPGo0Gn3/+eb/bXOiISnR0NFpaWnjVDxGRB3I6ncjLy4NKpYLZbIZGo+EEnzInSRJsNhvq6urgcDgwatSo8yZ1a21thdFovKzPb1kcUSkpKcHmzZuxbt26i26n1Wqh1WqHKRUREYlms9ngdDoRHR0NPz8/0XHoMvn6+sLHxwclJSWw2WzQ6XSDfi5ZrPWzatUqhIaG4sYbbxQdhYiIZOhS06yT/LjqZyb8J+90OrFq1SosWrQIarUsDvAQERGRTAhvBps3b0ZpaSnuu+8+0VGIiMhNVDR3oanDNmyvF+ivQaTJd9hej/5DeFGZO3cuZDKel4iI3EBFcxeufz4T3T3OYXtNnY8SW34zU2hZWb16NR555BE0NzcLyyCC8KJCREQ0EE0dNnT3OLF4VsKwFIeK5i68ui0fTR22Ab3ePffcg3feeee8++fNm4cNGzZc9LGxsbF45JFH8Mgjj/Ted+edd2LBggWX/fqDJbdCxKJCRERuKdLki7gQf9ExLuqGG27AqlWr+tw32KtXfX194evrfaefhA+mJSIi8lRarRYWi6XPLTAwEJIkYdmyZYiJiYFWq0VERAR++ctfAgBmzpyJkpIS/PrXv+5drBc4e6TDZDL1PveyZcswbtw4vP3224iJiUFAQAAeeughOBwO/N///R8sFgtCQ0Px9NNP98n0wgsvIC0tDf7+/oiOjsZDDz2E9vZ2AEBmZibuvfdetLS09L72smXLAJydy+zRRx9FZGQk/P39MWXKFGRmZg75PuQRFSIakIZ2K9bsKcGxihbk17bD6OuDCSMCcV1yKKYnmi/9BESETz/9FH/961+xdu1apKamorq6GkeOHAEArFu3Dunp6fj5z3+On/3sZxd9noKCAnz99dfYsGEDCgoK8P3vfx+FhYVITEzE9u3bsXv3btx3332YPXs2pkyZAuDsZcMvvfQS4uLiUFhYiIceegiPPfYYXnvtNVx99dVYsWIFnnjiCeTm5gIAAgICAAAPP/wwTp48ibVr1yIiIgKfffYZbrjhBhw7dgyjRo0asn3FokJEl6XH4cSqXUV4eUs+HJKExDA90iKNaO3uwcYT1Vi9uxizksxYdnMqRgTL+3A80XD54osvej/oz/n9738PnU4Hi8WC2bNnw8fHBzExMZg8eTIAICgoCCqVCnq9HhaL5aLP73Q68fbbb0Ov1yMlJQWzZs1Cbm4uvvrqKyiVSiQlJeG5557Dtm3beovKt8e9xMbG4k9/+hMeeOABvPbaa9BoNDAajVAoFH1eu7S0FKtWrUJpaSkiIiIAAI8++ig2bNiAVatW4ZlnnnHF7rogFhUiuqTuHgceei8bmbm1mD06DN+bEAWDzqf3+5Ik4WBJE/6xpxjz/pqF138yATOTQgUmJpKHWbNmYeXKlX3uCwoKQkdHB1asWIGRI0fihhtuwIIFC7Bw4cIBzycWGxsLvV7f+3VYWBhUKlWfydbCwsJQW1vb+/XmzZuxfPlynD59Gq2trbDb7eju7kZnZ2e/s/8eO3YMDocDiYmJfe63Wq0IDg4eUOaBYlEhoovqsjnw838cxL7CRvz3vGSMizadt41CocCk2CCkRRrxytZ83P/OQbx213jMTb34/waJPJ2/vz8SEhLOuz8oKAi5ubnYvHkzNm3ahIceegh//vOfsX37dvj4+FzgmS7su9sqFIoL3ud0nr2Uu7i4GDfddBMefPBBPP300wgKCsLOnTvx05/+FDabrd+i0t7eDpVKhUOHDkGlUvX53nePGLkaB9MSUb8kScIvPziM/UWN+O95SRcsKd+m81HhkdmjMDE2EA++m41tubUX3Z7Im/n6+mLhwoV46aWXkJmZiT179uDYsWMAAI1GA4fD4fLXPHToEJxOJ55//nlcddVVSExMRGVlZZ9tLvTaGRkZcDgcqK2tRUJCQp/bpU5PXSkeUSGifr2zuxibTtXgN3MTMSbSeFmPUauUeHjWKLxgz8Wv1h7GV7+8FlGBXEyOXK+iuUv2r2O1WlFdXd3nPrVajS+++AIOhwNTpkyBn58f3n33Xfj6+mLEiBEAzp7SycrKwn/9139Bq9UiJCTkiv4O5yQkJKCnpwcvv/wyFi5ciF27duH111/vs01sbCza29uxZcsWpKenw8/PD4mJibjrrrtw99134/nnn0dGRgbq6uqwZcsWjB07dkjX6mNRIaILOl7Rgqe/OoUbUi2YOCJoQI9VKRV4cEYC/mf9MTz0bjY+fnAqtGrVpR9IdBkC/TXQ+Sjx6rb8YXtNnY8Sgf6aAT9uw4YNCA8P73NfUlISnn32WTz77LNYsmQJHA4H0tLS8Pnnn/eO9/jDH/6AX/ziF4iPj4fVanXZDO7p6el44YUX8Nxzz2Hp0qWYPn06li9fjrvvvrt3m6uvvhoPPPAA7rzzTjQ0NODJJ5/EsmXLsGrVKvzpT3/Cb37zG1RUVCAkJARXXXUVbrrpJpdk649CcuP561tbW2E0GtHS0gKDwSA6DpHHsNmdmP9iFpwS8NTNqfBRDe4scUFdO576/AQWTY3F4zeluDgleYPu7m4UFRUhLi4OOp2u936u9SN//f3sgIF9fvOIChGdZ82eYhTVd+CZ29IGXVIAIN4cgDsmROPtXUW4NSPysk8fEV1KpMmXxcFLcDAtEfXR2GHDi5vzcF1yqEvmQ5mfZkFUoC/+57NjcDjd9gAuEQnCokJEffx10xk4JQl3TIh2yfOplUrcN20kjpS34IP9pS55TiLyHiwqRNSroK4d7+0rwW0ZUTD4Xv5cDpeSZNFjZpIZ/7fxNFq7e1z2vETk+VhUiKjX65kFMPlpMDc1zOXPfceEaHTbnPj7jiKXPzd5Pje+7sNruepnxqJCRADOXkWx7nAF5o+xXNEA2v4E+WswJyUMb+0oROMwXq1B7u3cLKudnZ2Ck9BAnfuZDWSm3QvhVT9EBAD4W1YhfH1UuD7Z9UdTzrl5XAS2nq7FG9sLsHTB6CF7HfIcKpUKJpOpd60aPz8/KBQKwanoYiRJQmdnJ2pra2Eymc6bcn+gWFSICA3tVqzdX4oFY8Phqxm6idkMOh/MH2PBO7uLcf+1I2HWa4fstchznJui/dsL65H8mUwml0yvz6JCRHh3bykkAPOGYRHB+Wnh+Pp4NdbsKcZv5iYN+euR+1MoFAgPD0doaCh6ejgY2x34+Phc8ZGUc1hUiLyc3eHEB/tLcXV8CAw6113p058ArRozksxYs6cED86Mh5+Gb0N0eVQqlcs+/Mh9cDAtkZfberoW1a3dmD06dNhec8EYC9q6e/DpofJhe00ick8sKkRe7r19pUgw+2OkOWDYXtOs12FyXBD+tqOIs9US0UWxqBB5sdKGTmSdqcP1o4fuSp/+3DQ2AqWNndh0smbYX5uI3AeLCpEX++BAKfw0KkyNDx721443ByAxLADv7i0Z9tcmIvfBokLkpRxOCZ8eKse0hBBo1WIGKF6fHIad+fUoru8Q8vpEJH8sKkReandBPWrbrLh2lFlYhqtGBiNAq+ZihUTULxYVIi/1WXYFIow6xJv9hWXQqJW4dlQIPjpYBqvdISwHEckXiwqRF+q02fH18WpMSwgRPh359aPD0NTZgw3Hq4XmICJ5YlEh8kIbT1Sjq8eBa0eFiI6CSJMvUsINeH8fT/8Q0flYVIi80LrsCoy26GHW60RHAQBMTzRjX1Ejypu4Qi4R9cWiQuRl6tut2JVfj2kJ4o+mnDM5NghatRLrD1eIjkJEMsOiQuRlNp44OxZkUlyQ4CT/4atRYXJcED45VA5J4ky1RPQfLCpEXuarY1VIjTAOywKEA3HtKDOKGzpxuKxZdBQikhEWFSIv0thhw96CRkyW0dGUc1LDDQj212BdNhcqJKL/YFEh8iIbT1RDgoRJsfIrKkqlAtMSQvCvnErY7E7RcYhIJlhUiLzIV8eqkBJugNFXXqd9zrk6Phit3XbszK8THYWIZIJFhchLNHXYsDu/QZanfc6JCfJDpMkXXxypEh2FiGSCRYXIS2w+VQOnJM/TPucoFApcNTIIG09Uo7uHU+oTEYsKkdfYfKoGo8ICYPLTiI5yUVNHhqDD5kBmLk//EJEMikpFRQV+/OMfIzg4GL6+vkhLS8PBgwdFxyLyKN09Dmw/U4cJMYGio1xSZKAvRgT74YujlaKjEJEMqEW+eFNTE6ZNm4ZZs2bh66+/htlsRl5eHgID5f9mSuROduXXo7vHiQkyPu3zbVfFBeOfRyrQabPDTyP0bYqIBBP6DvDcc88hOjoaq1at6r0vLi5OYCIiz7TpZA0iTDpEmnxFR7ksV40MxocHy5CZW4cFaeGi4xCRQEJP/fzrX//CxIkTcccddyA0NBQZGRn429/+1u/2VqsVra2tfW5EdHFOp4RNp2ow3g1O+5xjMeowItgPG49Xi45CRIIJLSqFhYVYuXIlRo0ahY0bN+LBBx/EL3/5S7zzzjsX3H758uUwGo29t+jo6GFOTOR+Dpc1o6Hdhokj3OO0zzkTRwRiy+laTv5G5OWEFhWn04nx48fjmWeeQUZGBn7+85/jZz/7GV5//fULbr906VK0tLT03srKyoY5MZH72XyqBkZfH4wKDRAdZUAmxQah3WrH7oJ60VGISCChRSU8PBwpKSl97hs9ejRKS0svuL1Wq4XBYOhzI6KL23qqFmOjjFAqFaKjDEhMkB/CDFpsPFEjOgoRCSS0qEybNg25ubl97jtz5gxGjBghKBGRZ6lo7kJuTRsyot1nfMo5CoUCk2KD8M2Jajickug4RCSI0KLy61//Gnv37sUzzzyD/Px8vP/++3jzzTexePFikbGIPMa207VQKoCxUUbRUQZlUmwQGjpsyC5tEh2FiAQRWlQmTZqEzz77DB988AHGjBmDP/7xj1ixYgXuuusukbGIPMa207VIsujhr3XPuUgSQgMQ6OeDb07w6h8ibyX83eumm27CTTfdJDoGkcfp7nFgV0E9bsuIEh1l0JQKBcZFB2LTyRr8z40pl34AEXkc4VPoE9HQ2FvYgO4eJzKiTaKjXJHxI0wobuhEYV276ChEJACLCpGH2na6Fma9FlGB7jEbbX/SIo3QqJTYcqpWdBQiEoBFhchDbcutQ3qUCQqFe12W/F1atQpjIg3YfIqXKRN5IxYVIg9UXN+B0sZOjHPz0z7nZMQE4mBxE1o6e0RHIaJhxqJC5IEyc2uhViqQGuEZkyKOjwmEQ5KQeYanf4i8DYsKkQfKPFOHZIseOh+V6CguEeSvwcgQf2w+ydM/RN6GRYXIw3T3OLC3oAFjo0yio7jUuGgTtp+p4yy1RF6GRYXIwxwobkS33Yl0Dxmfck56tAmt3XbklHGWWiJvwqJC5GG259YhyF+DaDe/LPm7EswB0OvUyMytEx2FiIYRiwqRh8k8U4exkUa3vyz5u5RKBdIijdh2mgNqibwJiwqRB6ls7kJ+bbvHjU85Jz3KhOOVrahrs4qOQkTDhEWFyIPsyKuDUnF2NldPdG4V6KwzPP1D5C1YVIg8SNaZesSbAxCgE77e6JAw+WkQb/ZHZi5P/xB5CxYVIg/hcErYkVfnsUdTzkmPMiErr56XKRN5CRYVIg9xrKIFrd12jx2fck56tAktXT04VtEiOgoRDQMWFSIPkXWmDn4aFeJD/UVHGVLx5gD4aVQcp0LkJVhUiDxE1pk6pEYYoFZ69j9r1b/XMGJRIfIOnv2ORuQl2rp7cLi0GWmRJtFRhkVapAmHS5vR1s3VlIk8HYsKkQfYU9AAhyT1Xr7r6dKjjHBIEvYUNIiOQkRDjEWFyAPsyKuHxaBDmEEnOsqwCDXoEG7UISuPp3+IPB2LCpEHyDpThzEeflnyd42JNCLrTL3oGEQ0xFhUiNxcWWMnSho7vea0zzljI40obexESUOH6ChENIRYVIjc3I68eigVQGqEQXSUYZUSYYBSAezM51EVIk/GokLk5rLy6jAqVA8/jWdOm98fP40ao0L12JnHokLkyVhUiNyYwylhV349xkR619GUc1IjDdhd0AAnp9Mn8lgsKkRu7Gh5M9q8YNr8/oyJMKKlqwcnq1pFRyGiIcKiQuTGdubVn5023xwgOooQo0IDoFUrOU6FyIOxqBC5say8s9Pmq5QK0VGEUKuUGB1u4DgVIg/GokLkptqt9n9Pm+9dlyV/15gIIw4UN6K7xyE6ChENARYVIje1t6ABdqfkteNTzhkTaYDV7kR2SZPoKEQ0BFhUiNzUjrw6hBm0XjNtfn+ig/xg9PXhOBUiD8WiQuSmtp+pw5gI7z7tAwBKhQIp4QbsKmBRIfJELCpEbqi8qRPFDZ1ef9rnnNQIA46Vt6C1u0d0FCJyMRYVIje000unze9PaoQRTgk4UNQoOgoRuRiLCpEbysqrQ3xoAPy13jVtfn/CDFqYA7TYXdAgOgoRuRiLCpGbOTttfgPSOD6ll0KhQEqEAbs4oJbI47CoELmZYxUtaOnq4fiU70iNMOB0dRsa2q2ioxCRC7GoELmZHWfqzk6bH+ovOoqspISfHa+zt5DjVIg8CYsKkZvZfubstPlqJf/5fltwgBYRJh128zJlIo/CdzoiN9LW3fPvafNNoqPIUko4x6kQeRoWFSI3srugAQ5JwtgoDqS9kJRwI4obOlHT2i06ChG5iNCismzZMigUij635ORkkZGIZG1HXh3CjTqvnza/Pyn/nldmDy9TJvIYwidhSE1NxebNm3u/VquFRyKSrawz9Rjj5aslX4zR1wfRgb7YU9CAWzMiRcchIhcQ3grUajUsFovoGESyV9LQgdLGTtwxIUp0FFkbHW7ggFoiDyJ8jEpeXh4iIiIwcuRI3HXXXSgtLe13W6vVitbW1j43Im+x/UwdVEoFUjnR20WlRhhR1tSFiuYu0VGIyAWEFpUpU6Zg9erV2LBhA1auXImioiJce+21aGtru+D2y5cvh9Fo7L1FR0cPc2IicTJz65Bs0cNXoxIdRdaSw/UAOE6FyFMILSrz58/HHXfcgbFjx2LevHn46quv0NzcjI8++uiC2y9duhQtLS29t7KysmFOTCSG1e7AnoIGjOX4lEsy6HwwItgPewtZVIg8gfAxKt9mMpmQmJiI/Pz8C35fq9VCq9UOcyoi8Q4VN6Grx4Gx0SbRUdwCx6kQeQ7hY1S+rb29HQUFBQgPDxcdhUhWtp+pQ6CfD0YE+YmO4hZSww2obO5GWWOn6ChEdIWEFpVHH30U27dvR3FxMXbv3o3bbrsNKpUKP/zhD0XGIpKdzNw6pEUaoVAoREdxC8nhBigA7OHpHyK3J7SolJeX44c//CGSkpLwgx/8AMHBwdi7dy/MZrPIWESyUtPajdyaNqTztM9lC9CqERviz3EqRB5A6BiVtWvXinx5IrewPbcOCoATvQ3QaIseewoaIEkSj0QRuTFZjVEhovNty61FQlgADDof0VHcSkqEEVUt3Shr5HwqRO6MRYVIxnocTuzIq0d6lEl0FLeTbNFDAfD0D5GbY1EhkrHskia0W+0Yx/EpA+avVSOO41SI3B6LCpGMbcutg9HXB3Eh/qKjuKXkcAN2F54dp0JE7olFhUjGtuXWYmyUEUoOBh2U1HADqlu6Ucr5VIjcFosKkUxVtXQht7oNGTztM2jJ4XooFVz3h8idsagQydT23DooFUBapEl0FLflpzk7n8q+okbRUYhokFhUiGRq6+lajArVI0AnqyW53M5oi6F3PhUicj8sKkQyZLU7sCOvHhkxJtFR3F5KuAHVrRynQuSuWFSIZGhvYSO6ehwYHxMoOorbS7KcHafCy5SJ3BOLCpEMbT1VA7Nei6hAX9FR3J5/77o/HKdC5I5YVIhkRpIkbD5Vi4xoE9eocZGUcI5TIXJXLCpEMpNX246K5i5k8LSPy4zmOBUit8WiQiQzW07VQqtWIiXcIDqKx0jmOBUit8WiQiQzW07VIC3SCI2a/zxdxU9zdt0fTvxG5H74TkgkI40dNmSXNvG0zxAYHW7AHq77Q+R2WFSIZGTLqRpIEjCe86e4XEq4ATWtVhQ3cJwKkTthUSGSkU0nazAqLAAmP43oKB6H86kQuScWFSKZ6O5xICuvDhN42mdI+GnUGGkO4DgVIjfDokIkE7vy69Hd48SE2CDRUTzWaIue41SI3AyLCpFMbDpZg3CjDhFGnegoHislwoC6NisK6ztERyGiy8SiQiQDTqeETadqMGFEIGejHUJJYQaOUyFyMywqRDJwuKwJDe02TBjB8SlDyVejQjzHqRC5FRYVIhnYeKIGJj8fJIbqRUfxeJxPhci9sKgQCSZJEr4+VoUJMYFQKnnaZ6ilRhjQ0G5Dfm276ChEdBlYVIgEO1nVirKmLkyO49U+wyExTA+1UoE9HKdC5BZYVIgE23i8Gv5aFRchHCY6HxUSQgOwO79edBQiugwsKkSCfX28GuOjA6FW8Z/jcEkJN2BPYSOcTo5TIZI7vjMSCVRQ14682nZM4mmfYZUaYUBLVw9OV7eJjkJEl8CiQiTQhuPV0KqVGBtlFB3FqySE6uGjUmB3AU//EMkdiwqRQF8erUJGjAlatUp0FK+iUSuRFKbngFoiN8CiQiRIcX0HTla14qq4YNFRvNLocAP2FTbC7nCKjkJEF8GiQiTIl8eqoFUrMS7GJDqKVxoTaUS71Y7jla2ioxDRRbCoEAnC0z5ijTT7w9dHxXEqRDLHokIkAE/7iKdWKpFs0WMX51MhkjUWFSIBeNpHHlIjjDhY3ITuHofoKETUDxYVIgE+P1LJ0z4ykBppgNXuRHZpk+goRNQPFhWiYZZf24bT1W24Oj5EdBSvFxPkB4NOjd35vEyZSK5YVIiG2b9yKuGvUWFctEl0FK+nVCiQEmHALg6oJZItFhWiYSRJEv6ZU4mJsUHw4do+spAaYcTRsha0dfeIjkJEF8B3SqJhdKyiBSWNnbg6nlf7yEVqhAEOScL+okbRUYjoAmRTVJ599lkoFAo88sgjoqMQDZl/5VTC6OuD1Aiu7SMXFoMO5gAtdvIyZSJZkkVROXDgAN544w2MHTtWdBSiIeN0Svj8aCUmxwVBpVSIjkP/plAoMCbSgB15LCpEciS8qLS3t+Ouu+7C3/72NwQGBoqOQzRk9hY1oKbVimm82kd2xkQakV/bjprWbtFRiOg7BlVUCgsLXRZg8eLFuPHGGzF79uxLbmu1WtHa2trnRuQu1h+uQJhBi8SwANFR6DvG/PtU3E4eVSGSnUEVlYSEBMyaNQvvvvsuursH/z+QtWvXIjs7G8uXL7+s7ZcvXw6j0dh7i46OHvRrEw2n7h4HvjpWjWkJIVAoeNpHbgy+PogL8ec4FSIZGlRRyc7OxtixY7FkyRJYLBb84he/wP79+wf0HGVlZfjVr36F9957Dzqd7rIes3TpUrS0tPTeysrKBhOfaNhtPlWDdqsd1/C0j2ylRhiwI68OkiSJjkJE3zKoojJu3Di8+OKLqKysxNtvv42qqipcc801GDNmDF544QXU1dVd8jkOHTqE2tpajB8/Hmq1Gmq1Gtu3b8dLL70EtVoNh+P8tTe0Wi0MBkOfG5E7+Cy7AgmhAQg3+YqOQv1IizSivt2GMzXtoqMQ0bdc0WBatVqN22+/HR9//DGee+455Ofn49FHH0V0dDTuvvtuVFVV9fvY66+/HseOHUNOTk7vbeLEibjrrruQk5MDlYproJBnaOywYfuZOg6ilblkiwE+KgV25F36P1pENHyuqKgcPHgQDz30EMLDw/HCCy/g0UcfRUFBATZt2oTKykrccsst/T5Wr9djzJgxfW7+/v4IDg7GmDFjriQWkaz8K6cCEsBJ3mROo1Yi2cLLlInkRj2YB73wwgtYtWoVcnNzsWDBAqxZswYLFiyAUnm298TFxWH16tWIjY11ZVYit/TJoXJkRJtg8PURHYUuYWyUEZ8eKkd3jwM6Hx7VJZKDQRWVlStX4r777sM999yD8PDwC24TGhqKt956a0DPm5mZOZg4RLKVW92G45WtWDInUXQUugxjo0x4b18pDhQ34tpRZtFxiAiDLCqbNm1CTExM7xGUcyRJQllZGWJiYqDRaLBo0SKXhCRyV59ml0OvUyODKyW7hehAXwT5a5B1po5FhUgmBjVGJT4+HvX155/HbWxsRFxc3BWHIvIEdocTn2aXY1p8CNRcKdktKBQKpEUasf0MB9QSycWg3j37m2egvb39sudEIfJ0O/Lq0dBuw/RE/s/cnaRHGXGmph1VLV2ioxARBnjqZ8mSJQDO/q/jiSeegJ+fX+/3HA4H9u3bh3Hjxrk0IJG7+vBgGWKC/BAb7HfpjUk2xkQaoQCw40w9fjCJs18TiTagonL48GEAZ4+oHDt2DBqNpvd7Go0G6enpePTRR12bkMgNNbRbsflkDX40JYZT5rsZvc4H8aEB2H6mjkWFSAYGVFS2bdsGALj33nvx4osvcmZYon58drgCCgVwTQIneXNHY6OM2HyyBnaHk+OLiAQb1L/AVatWsaQQ9UOSJHx4oAwTRgRCr+PcKe5oXJQJrd125JQ1i45C5PUu+4jK7bffjtWrV8NgMOD222+/6Lbr1q274mBE7upIeQvyatvxvfHJoqPQIMWbA2DQqbEttxYTY4NExyHyapddVIxGY++5dqPROGSBiNzdhwfKEBKgQVok/524K6VSgbFRJmw9XYv/nsfCSSTSZReVVatWXfDPRPQfHVY7/pVTgXmpFiiVHETrzsZFm/DKtnxUt3TDYuS0C0SiDGqMSldXFzo7O3u/LikpwYoVK/DNN9+4LBiRO/riaCU6bQ7MTAoVHYWuUHqUCUoFkJlbKzoKkVcbVFG55ZZbsGbNGgBAc3MzJk+ejOeffx633HILVq5c6dKARO7k/f2lSI82wqzXio5CVyhAp8aoMD22sagQCTWoopKdnY1rr70WAPDJJ5/AYrGgpKQEa9aswUsvveTSgETu4lRVK46UteC6pDDRUchFxkWZsCOvHja7U3QUIq81qKLS2dkJvV4PAPjmm29w++23Q6lU4qqrrkJJSYlLAxK5i7X7S2Hy9UHGCJPoKOQiGTEmdNoc2F/UKDoKkdcaVFFJSEjA+vXrUVZWho0bN2Lu3LkAgNraWs6vQl6py+bAusMVmJ5ohlrJCcI8RUyQH0ICNNh8qkZ0FCKvNah31CeeeAKPPvooYmNjMWXKFEydOhXA2aMrGRkZLg1I5A6+OFqJtm47rkvmIFpPolAoMD4mEJtO1vS7GCsRDa0BTaF/zve//31cc801qKqqQnp6eu/9119/PW677TaXhSNyF+/uK0F6lBFhBl7G6mkmjAjENydrkFvThmQLjxgTDbdBFRUAsFgssFgsfe6bPHnyFQcicjcnKltwpKwFS+Ykio5CQ2B0uAG+PipsPlnDokIkwKCKSkdHB5599lls2bIFtbW1cDr7jogvLCx0STgid/D+vlIE+WswPiZQdBQaAj4qJcZGGfHNyRo8fN0o0XGIvM6gisr999+P7du34yc/+QnCw8O5jD15rXarHZ8drsANYyxQcSZajzVhRCBeyyxAbWs3Qnl6j2hYDaqofP311/jyyy8xbdo0V+chcivrD1egu8eB6zgTrUcbF312ltrNp2rxoykxouMQeZVBXfUTGBiIoCCuKEreTZIk/GNPCSaMCERwAGei9WR6nQ+SLQZsPFEtOgqR1xlUUfnjH/+IJ554os96P0Te5mBJE3Jr2jB7NGei9QaTYgOxK78erd09oqMQeZVBnfp5/vnnUVBQgLCwMMTGxsLHx6fP97Ozs10SjkjO/rGnGOFGHcZEGkVHoWEwKTYI7+wpwbbTtbhlXKToOEReY1BF5dZbb3VxDCL3UtdmxVfHqvHDyTFQcjC5VwgO0CLe7I8Nx6tZVIiG0aCKypNPPunqHERu5aODZVApFZieaBYdhYbRxNgg/DPn7ABqnY9KdBwirzDoRUmam5vx97//HUuXLkVj49kFu7Kzs1FRUeGycERyZHc48e7eEkwdGYwA7aDnTCQ3NDk2CN09TmSdqRMdhchrDKqoHD16FImJiXjuuefwl7/8Bc3NzQCAdevWYenSpa7MRyQ7m0/VoqqlG3NTLZfemDxKhMkXUYG+2MCrf4iGzaCKypIlS3DPPfcgLy8POt1/Jj9asGABsrKyXBaOSI7e2V2MpDA94kL8RUchASbHBmHTiRrY7M5Lb0xEV2xQReXAgQP4xS9+cd79kZGRqK7m/zTIc+XVtGFPYQPmpPCSZG911chgtFnt2JnP0z9Ew2FQRUWr1aK1tfW8+8+cOQOzmYMLyXOt2VMCk68PpsRxwkNvFRV49vTPF0erREch8gqDKio333wz/vCHP6Cn5+zERwqFAqWlpfjtb3+L733vey4NSCQXrd09+ORQOa4bHQq1atDj0MnNKRQKTIkLwjcnamC1O0THIfJ4g3q3ff7559He3g6z2Yyuri7MmDEDCQkJ0Ov1ePrpp12dkUgWPj5YDpvDyZloCVPigtFutWPHmXrRUYg83qCurTQajdi0aRN27dqFI0eOoL29HePHj8fs2bNdnY9IFpxOCe/sLsaUuCAE+mlExyHBooP8EB3oi8+PVmI2xysRDakBFxWn04nVq1dj3bp1KC4uhkKhQFxcHCwWCyRJgoKzdJIHyjxTi9LGTtx/TZzoKCQTU0YG46tjVZz8jWiIDejUjyRJuPnmm3H//fejoqICaWlpSE1NRUlJCe655x7cdtttQ5WTSKhVu4oRb/ZHQmiA6CgkE1fHB6PT5sDW07WioxB5tAEdUVm9ejWysrKwZcsWzJo1q8/3tm7diltvvRVr1qzB3Xff7dKQRCLl17ZjR149HpwRzyOG1Cvc6IsEsz/WH67AgrRw0XGIPNaAjqh88MEH+P3vf39eSQGA6667Dr/73e/w3nvvuSwckRys3l0Ek68PpsYHi45CMjM1PgTbcmvR0tkjOgqRxxpQUTl69ChuuOGGfr8/f/58HDly5IpDEclFS+fZS5KvHx0GH16STN8xNT4YDqeEr49zThWioTKgd97GxkaEhfU/wj0sLAxNTU1XHIpILtYeKIXDKWH26FDRUUiGAv00SI0wYn1OpegoRB5rQEXF4XBAre5/WItKpYLdbr/iUERyYHc48c7uYkwdGQwTL0mmfkxLCMa+wgZUtXSJjkLkkQY0mFaSJNxzzz3QarUX/L7Vah3Qi69cuRIrV65EcXExACA1NRVPPPEE5s+fP6DnIRoKm07WoLKlGw9fN0p0FJKxSbFBeHtnMdYfrsSDM+NFxyHyOAMqKosWLbrkNgO54icqKgrPPvssRo0aBUmS8M477+CWW27B4cOHkZqaOpBoRC73951FGB3OVZLp4vw0akyMDcQnh8rwwIyRvDKMyMUUkiRJokN8W1BQEP785z/jpz/96SW3bW1thdFoREtLCwwGwzCkI29xpKwZt7y6C0vmJGJSLBcgpIs7UtaMZzecxj8XT0N6tEl0HCLZG8jn96Cm0B8KDocDH3/8MTo6OjB16tQLbmO1WvucXrrQCs5ErvDWzkKEGbSYEBMoOgq5gbRII4L8Nfg0u5xFhcjFhF9veezYMQQEBECr1eKBBx7AZ599hpSUlAtuu3z5chiNxt5bdHT0MKclb1DZ3IUvj1XjhlQLlEoexqdLUyoVmBYfjH/mVHJFZSIXE15UkpKSkJOTg3379uHBBx/EokWLcPLkyQtuu3TpUrS0tPTeysrKhjkteYN39hRDq1ZiRiIvSabLNz3RjJauHmzjlPpELiX81I9Go0FCQgIAYMKECThw4ABefPFFvPHGG+dtq9Vq+73iiMgVOqx2vL+vFLOSQuGr4UJzdPmiAv2QYPbHhwfKcMMYTqlP5CrCj6h8l9PpHPBlzkSu8tHBMnRY7bhhjEV0FHJD0xNDsf1MHapbukVHIfIYQovK0qVLkZWVheLiYhw7dgxLly5FZmYm7rrrLpGxyEs5nBLe2lmEq0YGIySAR+5o4KYlBMNHpcSn2eWioxB5DKFFpba2FnfffTeSkpJw/fXX48CBA9i4cSPmzJkjMhZ5qW9OVKO8qYsr4dKg+WnUmBwXhA8PlMHplNXMD0RuS+gYlbfeekvkyxP18eaOQowO1yPeHCA6CrmxWUmh+MMXJ7GvqJErbhO5gOzGqBCJcKikEYdLm3k0ha5YskWPcKMOaw+Uio5C5BFYVIgAvJlViEiTL8Zzgje6QgqFAjMTzfj6WDWaO22i4xC5PRYV8npF9R345kQNFqSFQ8l1WsgFZiSFwilJ+DS7QnQUIrfHokJe762dhTD4+uCahBDRUchDGH19MDE2EO/tK4HMllMjcjssKuTVGtqt+PhgOeamhEGj5j8Hcp3rk8NQWNeB/UWNoqMQuTW+M5NXe2dPCQBgTkqY4CTkaVIjDAg36vD+Pg6qJboSLCrktTptdryzuxizkkKh1/mIjkMeRqFQ4LrkUHx1vAqNHRxUSzRYLCrktT46UIa27h4sSON0+TQ0pieaAQAfH+QCqkSDxaJCXsnucOLNHYW4Oj4EZr1OdBzyUAadD6aODMY/9pZwplqiQWJRIa/05bEqVDZ346axnOCNhtaclDCUN3Vhe16d6ChEbolFhbyOJEl4LbMA6VFGjAj2Fx2HPFy8OQBxIf74x78HbhPRwLCokNfJPFOH3Oo23JweIToKeQGFQoHZo8Ow7XQtyho7RcchcjssKuR1Vm4rwKjQAIwON4iOQl7i6vhg+GlUeI+XKhMNGIsKeZVDJY3YX9yIhekRUHC6fBomOh8VZiSasfZAKbp7HKLjELkVFhXyKiszCxBp8sWEEVx8kIbXnBQLmjt78PmRStFRiNwKiwp5jdzqNmw+VYuF6Vx8kIafxajDuGgT3tldzPV/iAaARYW8xsrMfIQEaDAtnosPkhhzU8JwvLIVOWXNoqMQuQ0WFfIKZY2d+PxIFW5MC4daxV97EiM92oQwgxardxeLjkLkNviOTV7hjawC+GtVmJUcKjoKeTGlQoE5oy348mgVatu6RcchcgssKuTxatu68dGBcsxLtUCrVomOQ15uZpIZKqUCH+zj+j9El4NFhTzeWzuKoFYpMC+Viw+SeP5aNa5JCMG7+0pgsztFxyGSPRYV8mjNnTb8Y28J5qSEwV+rFh2HCAAwL9WCujYrNpyoFh2FSPZYVMijrd5dDIdTwvwxXHyQ5CM6yA+pEQas3lUkOgqR7LGokMdqt9rx9q4izEoKhdHXR3Qcoj7mpViQXdqM4xUtoqMQyRqLCnms9/aWoNPqwE1jeTSF5Gf8iECYAzRYvZtHVYguhkWFPFJ3jwNv7ijE9EQzggO0ouMQnUelPLuq8r9yqtDQbhUdh0i2WFTII314oAxNHTbcnB4hOgpRv87N67P2AC9VJuoPiwp5HJvdiZWZBbg6PgRhBp3oOET90ut8cHV8MP6xpwR2By9VJroQFhXyOJ8dLkd1azduGcejKSR/88ZYUN3ajU0na0RHIZIlFhXyKHaHE69uK8Dk2CBEBfqJjkN0SbHB/ki26Ln+D1E/WFTIo3xxtAqljZ24NSNSdBSiyzY3JQz7ihpxurpVdBQi2WFRIY/hdEp4eWsexseYEBfiLzoO0WWbFBeEQD8frNlTIjoKkeywqJDH2HCiGgV1Hbh1HI+mkHtRK5W4fnQY1mWXo6WzR3QcIllhUSGPIElnj6akRRoxKkwvOg7RgF2fHAq7Q8LHh3ipMtG3saiQR9h6uhanqto4NoXclslPg8lxQfjHnhI4nZLoOESywaJCbk+SJLy0NQ/JFj1GW3g0hdzX3BQLSho7sSO/XnQUItlgUSG3tyu/AUfKWnDruEgoFArRcYgGLTEsALHBfljDS5WJerGokNt7eWse4s3+GBtlFB2F6IooFArMTgnD1tO1KGvsFB2HSBZYVMit7S9qxL6iRtzCoynkIabFh8BPo8J7+0pFRyGSBRYVcmuvbM1DTJAfJowIFB2FyCV0Pipcm2jGhwdKYbU7RMchEo5FhdzWkbJmZOXV49ZxEVDyaAp5kNmjw9DU2YMNx6tFRyESTmhRWb58OSZNmgS9Xo/Q0FDceuutyM3NFRmJ3MgrW/MRYdRhSlyw6ChELhVp8kVqhIEz1RJBcFHZvn07Fi9ejL1792LTpk3o6enB3Llz0dHRITIWuYHT1a3YdKoGN4+LhFLJoynkeeaMDsOhkiau/0NeTy3yxTds2NDn69WrVyM0NBSHDh3C9OnTBaUid/Dq1nyY9VpMS+DRFPJME2IDEejng3f3luBPt6aJjkMkjKzGqLS0tAAAgoKCLvh9q9WK1tbWPjfyPoV17fjiaBUWjg2HWimrX2Eil1ErlZiVFIp12RVot9pFxyESRjbv8k6nE4888gimTZuGMWPGXHCb5cuXw2g09t6io6OHOSXJwcrMApj8fDAjMVR0FKIhNSs5FN09Dvwrp1J0FCJhZFNUFi9ejOPHj2Pt2rX9brN06VK0tLT03srKuHiXtylv6sS6wxW4MS0CGrVsfn2JhkRIgBbjok14bx8H1ZL3ksU7/cMPP4wvvvgC27ZtQ1RUVL/babVaGAyGPjfyLm9mFcJPo8L1o3k0hbzD9clhOFHZiqPlzaKjEAkhtKhIkoSHH34Yn332GbZu3Yq4uDiRcUjmatu6sXZ/GW5ItUDnoxIdh2hYjIs2ISRAg/c5Uy15KaFFZfHixXj33Xfx/vvvQ6/Xo7q6GtXV1ejq6hIZi2Tq7Z3FUCkVmJtqER2FaNgolQrMSgrF+pwKtHb3iI5DNOyEFpWVK1eipaUFM2fORHh4eO/tww8/FBmLZKilswf/2FOMOSlhCNAKvaqeaNjNTAqFze7koFrySkLf8SVJEvny5EbW7ClGj0PC/DE8mkLeJ8hfg/ExgXh/Xyl+fNUI0XGIhpUsBtMSXUynzY63dhVhRpIZJj+N6DhEQsxKDsXJqlYcK28RHYVoWLGokOyt3V+G1q4eLBwbLjoKkTDjokwIDtDg/f0cVEvehUWFZM1md+LNrEJMSwiBWa8THYdIGKVSgRmJZvwzpwIdnKmWvAiLCsna+pwKVLd2Y+HYCNFRiISblRSKLpsDnx/hoFryHiwqJFtOp4SVmQWYOCIQ0UF+ouMQCRcSoEV6tBFrD3BWbvIeLCokW9+crEZRfQduTufRFKJzZiaFIqesGbnVbaKjEA0LFhWSJUmS8Nq2AqRGGDAqTC86DpFsTIgJhNHXBx/yqAp5CRYVkqU9BQ04WtHCsSlE36FWKXFNQgjWHS6H1e4QHYdoyLGokCy9llmA2GA/jI0yio5CJDuzkkLR3NmDTSdrREchGnIsKiQ7xytasDO/HgvTI6BQKETHIZKdyEBfJFn0WLufp3/I87GokOyszCxAmEGLKXHBoqMQydbMRDN25dejvKlTdBSiIcWiQrJS0tCBr49XYUFaOFRKHk0h6s9VI4Oh81Hhk0PloqMQDSkWFZKVN7MKodf5YGZiqOgoRLKm81HhqpFB+OhgGZxOLvBKnotFhWSjrs2Kjw+WY25KGDRq/moSXcrMpFBUNndjd0GD6ChEQ4afBiQbq3cXQakE5qZYREchcgujQgMQafLFhwe4UCF5LhYVkoV2qx1r9pTguqRQBOjUouMQuQWFQoGZSWZsPFGD5k6b6DhEQ4JFhWRh7f5SdNocWJAWLjoKkVu5JiEEdqcT/+JCheShWFRIOJvdib/vKMK0+GAEB2hFxyFyKyY/DcbHBHJKffJYLCok3OdHKlHd2o2bOF0+0aDMTArFicpWnKhsER2FyOVYVEgop1PC69sLMD7GhOggP9FxiNzSuGgTAv188PFBzqlCnodFhYTKPFOLvNp2Hk0hugIqpaJ3ocLuHi5USJ6FRYWEej2zEKPCApBs0YuOQuTWZiSForXLzoUKyeOwqJAw2aVN2F/ciIVpXHyQ6EpFmnyRFKbHxwc5qJY8C4sKCfPG9gJEGHWYMCJQdBQijzAj0YwdefWoaO4SHYXIZVhUSIiCunZ8c6IGC8aGQ8nFB4lc4qqRwdD6KPEpFyokD8KiQkL8LasQJj8fXJtgFh2FyGP4alSYEhfMhQrJo7Co0LCrbe3Gp9nluCHVwsUHiVxsZpIZ5U1d2FvIhQrJM/BTgobdW7uK4KNSYnZKmOgoRB4nKUyPCKOOM9WSx2BRoWHV2t2D9/aW4vrkUPhpuPggkaspFArMSArF18er0dLZIzoO0RVjUaFh9d7eUljtDszn4oNEQ2b6qBA4JAnrcypERyG6YiwqNGy6exz4+85CXDvKjEA/jeg4RB7r7EKFJnywvxSSxEG15N5YVGjYfHKoHE0dNizkdPlEQ25mUihOV7fheEWr6ChEV4RFhYaF3eHEG9sLMDkuCBajTnQcIo+XHmVCkL8Gaw+Uio5CdEVYVGhYfHW8GmVNXbg5PVJ0FCKvoFIqMH2UGf/MqUSnzS46DtGgsajQkJMkCa9ty8fYSCPiQvxFxyHyGrOSzOiw2vHF0SrRUYgGjUWFhtzW07U4Xd2GWzJ4NIVoOIUadEiLMuKD/Tz9Q+6LRYWGlCRJeHlrPpLC9Bht0YuOQ+R1rksOxeHSZuRWt4mOQjQoLCo0pPYWNiKnrBm3jIuAQsHFB4mG24QRgTD5+fCoCrktFhUaUq9uy0dssB/GRZtERyHySmqlEtNHmfFpdjm6exyi4xANGIsKDZnDpU3YmV+Pm9MjeTSFSKDrkkPR1s1BteSehBaVrKwsLFy4EBERZ08LrF+/XmQccrGXtuQhKtAXU+KCREch8mphBh3So4z4x95i0VGIBkxoUeno6EB6ejpeffVVkTFoCByvaMG23DrcMi4SSiWPphCJNnt0GI6UteB4RYvoKEQDInT52vnz52P+/PkiI9AQeWlLHixGHaaODBYdhYgAZMQEIjhAg/f2lWD57WNFxyG6bG41RsVqtaK1tbXPjeTnVFUrvjlZg1vHRUDFoylEsqBSKnBdUig+O1yB1u4e0XGILptbFZXly5fDaDT23qKjo0VHogt4cfMZhBm0mJYQIjoKEX3LrORQ9DgkfHqoXHQUosvmVkVl6dKlaGlp6b2VlZWJjkTfcaKyBRtO1OC2jEiolW7160Xk8QL9NJgSF4TVu4vhdEqi4xBdFrf6JNFqtTAYDH1uJC8rNp0dm3JNgll0FCK6gHmpFpQ0dGJ7Xp3oKESXxa2KCsnbsfIWbDpVg9szIjk2hUimRoUGIN7sj9W7ikVHIbosQotKe3s7cnJykJOTAwAoKipCTk4OSks51bM7+ss3uYgw+eLqeI5NIZIrhUKBuSkWbD9Th4K6dtFxiC5JaFE5ePAgMjIykJGRAQBYsmQJMjIy8MQTT4iMRYOwv6gR28/U4Y4JUTyaQiRzU+ODYfT14VEVcgtC51GZOXMmJIkDutydJEl4bsNpjAzxx2TOQkskez4qJWaPDsPHB8uwZE4iAv01oiMR9YtjVOiKZebW4VBJE+6YGA0l1/QhcgtzU8PglIB395aIjkJ0USwqdEUczrNHU0aH65EeZRQdh4guk0Hng+mJZqzeXcxVlUnWWFToiqzLLsfp6jb8aHIMV0gmcjM3poWjscOGzw5XiI5C1C8WFRq0LpsDf/kmF1eNDEJCqF50HCIaIItRh0mxQXgjqwAOTgBHMsWiQoP29q4iNLTb8F+TYkRHIaJBunlcBIrrO/HlsSrRUYguiEWFBqW2rRuvbcvH7JQwhBl0ouMQ0SDFmwMwLtqIl7fkcVp9kiUWFRqUP2/IhVKhwPcyokRHIaIrdFtGFPJq2/HNyWrRUYjOw6JCA3akrBkfHyrH9ydGIUAndCoeInKBxDA9UiMMeHFLHue2ItlhUaEBcTolPPmvE4gJ8sP1yWGi4xCRi9yeEYlTVW3YcJxHVUheWFRoQD7NLkdOWTMWTR3BqfKJPEhKhBFjo4z4v425sDucouMQ9WJRocvW1GHD01+dwjUJIUiJ4ORuRJ7mvybFoKi+A58cKhcdhagXiwpdtuVfn0KP3Ym7pvByZCJPFBfij6vjg/HXzWc4Wy3JBosKXZYDxY346GA57pwUA5MfFzAj8lR3TIhGfbsNf99RKDoKEQAWFboM3T0OPPbJUYwKC8D1o0NFxyGiIWQx6nBDqgWvbMtHZXOX6DhELCp0aX/ddAblTZ34xbXxXB2ZyAvcPj4Svj4qPP3lKdFRiFhU6OIOlzbhbzsK8f3xUYgM9BUdh4iGgZ9GjR9OjsGXx6qwO79edBzyciwq1K8umwO/+fgI4kL8cePYCNFxiGgYXZMQgiSLHr//7Bi6bBxYS+KwqFC/nv7qJCqauvDgjATOmULkZRQKBX5+7UhUNnfjzxtzRcchL8aiQhe06WQN3t1birumjOApHyIvFWHyxQ8mRmPVriLsL2oUHYe8FIsKnaeqpQuPfXIEE0cEYjav8iHyavPHWJBo0eM3H+WgpatHdBzyQiwq1IfN7sRD72VDqVDgZ9NHQsGrfIi8mlKpwIMz4tHYacNvPsqB08lFC2l4sahQH8u/PoVj5S341fWjYND5iI5DRDIQZtDhoRkJ2HyqFq9nFYiOQ16GRYV6rT9cgVW7ivHjq0ZgVJhedBwikpHxIwJx67hI/GVjLjadrBEdh7wIiwoBAA6VNOGxT45iemII5qaEiY5DRDJ0x4QoTIwNwuL3srGvsEF0HPISLCqE8qZO/HzNQcSb/XH/NRyXQkQXplQq8PCsBCSGBeCn7xzEsfIW0ZHIC7CoeLnGDhvufns/1CoFHpmTCB8VfyWIqH8+KiWWzElCuFGHO9/cgx15daIjkYfjp5IXa7facc/b+9HQbsNvb0jm4Fkiuiy+GhV+v2A0EsP0uGfVAXx0sEx0JPJgLCpeqsvmwM/eOYj8unb89oZkhBs5qRsRXT6djwq/mZuIGYlmPPbJUfzqg8OcZ4WGBIuKF+q02XHf6gPILm3Co3OTEBfiLzoSEbkhtVKJn107Eg/PSsDmUzWY99csfHKoHA7OtUIupJAkyW1/o1pbW2E0GtHS0gKDwSA6jltot9rx09UHcKS8Gb+dl4zkcO43IrpydW1WvLevBPuKGpEYFoB7p8VhYXoEArRql71Gl82BujYr6jusaOqwod1qR4fVgR6HEw6nBKUC0KhV8NUoYfLVwOTnA4tRh1C9juuVycxAPr9ZVLxIbVs37l11AEX1HfjveUlItnCfEZFr5de2YV12BY6UN0OjVmL6KDOmxgcjIyYQscF+MPlpLvg4m92Jhg4rqlq6Ud3SjfKmTlQ0daG8uQtljZ2obulGa7f9vMcpAKhVCigVCkgS0ONw4rsfamqlAhEmXySEBiAxTI8xkQakRRoRE+THqxwFYVGh8xTUtWPR2/vRYbXjtzckY0QwT/cQ0dBpaLdiR149jlY0I6+2HXbH2Y8aP40KAVo1dD4qSJIEm8OJtm47Om2OPo/39VEhRK9BiL8WwQFaBAdoEOyvgclPA4NODb3OB34aFbRqZZ+yce45O6wOtHX3oKHDhoZ2K6pbrahs7kJZUyca2m0AgCB/DSaMCMRVI4NxdXwwksL0UPLIy7BgUaE+tpyqwa/W5sDk54PH5iXDrNeKjkREXsRmd6KiuQs1rd2ob7eiu8cBq90JpUIBtUoBXx8V9DofGHRqBPlrEOyvhb9WNWRHO1q7elBQ14682nacqWnDmZo29DgkBAdoMCPRjJlJoZiRaIbRl1dCDhUWFQIAOJwSXt6ahxc352HCiEA8ODMefhrXnS8mIvIENrsTZ2racLS8GUfKW1Da2Am1UoHJcUGYl2rB3NQwXhnpYiwqhIrmLjyy9jAOlTTh9vFRuC0jEkqeiyUiuqT6diuyS5twqKQJJypb4XBKGBdtwo1p4ZifZkFUoJ/oiG6PRcWLSZKEDw+U4ekvT0Hro8TimQm8soeIaJA6rHZklzZhf1EjjpQ3o8chIT3aiIVjIzA/LRyRJh5pGQwWFS+VV9OGJ/55AnsKGzAj0YwfXzXCpZcGEhF5sy6bA9mlTdhb2NBbWjJizh5pWZAWjgiWlsvGouJlGjtseGlLHv6xpwQheg3umxaHsVEm0bGIiDxWp82O7NJm7C1owNGKvqWFR1oujUXFSzR12PD3nYVYtasYkgTcmhGJ+WMsXFiQiGgYddrsOFTS9/RQWqQRN4yxYG5KGBJCAzhfy3ewqHi4wrp2rNpVjI8PlUGSgHmpFtw4NpyLChIRCdZpsyOnrBkHihuRU9aM7h4nRgT5YXZKGK5LDsWk2CBo1PzPJIuKB+qw2vHNyWqs3V+GfUWNMPr6YPboMMxNCYOB1/oTEcmOze7E8coWHCppQk5ZMxo7bPD1UeHqhGBMH2XG1fHBXnu0xe2Kyquvvoo///nPqK6uRnp6Ol5++WVMnjz5ko/z9KLS0G5FVl4dNp6owbbTtbDanUiNMGBGohlT4oLZyomI3IQkSShp7MSRsmYcq2hBbnUb7E4JIQEaTIkLxuS4IEwYEYhkix5qLzh971ZF5cMPP8Tdd9+N119/HVOmTMGKFSvw8ccfIzc3F6GhoRd9rKcVlZrW7rOHDIsasbewAScqWyEBSDD7Y3JcMK4aGQSzXic6JhERXaHuHgfO1LThRGUrcqvbUFDXDrtTgs5HidQII9IijRgTaUSyRY+E0ADofFSiI7uUWxWVKVOmYNKkSXjllVcAAE6nE9HR0fh//+//4Xe/+91FH+uORcXplFDXbkVJQyeK6ztQUNeO3Oo2nKxqRW2bFQBgDtAi2aJHaqQB6VGmfhfxIiIiz2CzO1FU34G82jYU1neguL4DVS3dAAClAogK9ENCaABig/0xItgP0UG+iDT5IcKkQ4BW7Xanjwby+S10kg2bzYZDhw5h6dKlvfcplUrMnj0be/bsOW97q9UKq9Xa+3VLSwuAs3/hoSJJEpwSYHeeXUa8xyHB4XDC9u9bj92Jrh4Hum1OdNod6LTa0d5tR5u1B61ddjR19qC504baNivq/n2zO//TDc0BGkSYfDE+3BdxY4IRF+KHIP9vrcXjsKKjzXqBZERE5Emi/IGoOD1mxekBAF09DlQ0d6K8uQvVzd2oqm/C8eIq1LfZ+nyOaH2UMAdoEeyvQVCABoG+Ghj9fKDXqaHXqeGvUcNXo4a/9uwijjq1CjofJXzUSmhUSviolFCrlFArz669pFIooFIqhrT8nPvcvpxjJUKLSn19PRwOB8LCwvrcHxYWhtOnT5+3/fLly/HUU0+dd390dPSQZRxqZQCyRYcgIiK3li86wCC1tbXBaDRedBu3mrZ06dKlWLJkSe/XTqcTjY2NCA4ORltbG6Kjo1FWVuY2p4E8SWtrK/e/INz34nDfi8N9L9aV7n9JktDW1oaIiIhLbiu0qISEhEClUqGmpqbP/TU1NbBYLOdtr9VqodVq+9xnMpkAoPcQlcFg4C+tQNz/4nDfi8N9Lw73vVhXsv8vdSTlHKHXQGk0GkyYMAFbtmzpvc/pdGLLli2YOnWqwGREREQkB8JP/SxZsgSLFi3CxIkTMXnyZKxYsQIdHR249957RUcjIiIiwYQXlTvvvBN1dXV44oknUF1djXHjxmHDhg3nDbC9FK1WiyeffPK8U0M0PLj/xeG+F4f7Xhzue7GGc/8Ln0eFiIiIqD+eP08vERERuS0WFSIiIpItFhUiIiKSLRYVIiIiki23KiqvvvoqYmNjodPpMGXKFOzfv/+i269YsQJJSUnw9fVFdHQ0fv3rX6O7u3uY0nqWgez7np4e/OEPf0B8fDx0Oh3S09OxYcOGYUzrObKysrBw4UJERERAoVBg/fr1l3xMZmYmxo8fD61Wi4SEBKxevXrIc3qqge7/qqoq/OhHP0JiYiKUSiUeeeSRYcnpiQa679etW4c5c+bAbDbDYDBg6tSp2Lhx4/CE9TAD3fc7d+7EtGnTEBwcDF9fXyQnJ+Ovf/2ry/K4TVH58MMPsWTJEjz55JPIzs5Geno65s2bh9ra2gtu//777+N3v/sdnnzySZw6dQpvvfUWPvzwQ/z+978f5uTub6D7/vHHH8cbb7yBl19+GSdPnsQDDzyA2267DYcPHx7m5O6vo6MD6enpePXVVy9r+6KiItx4442YNWsWcnJy8Mgjj+D+++/nG/YgDXT/W61WmM1mPP7440hPTx/idJ5toPs+KysLc+bMwVdffYVDhw5h1qxZWLhwId93BmGg+97f3x8PP/wwsrKycOrUKTz++ON4/PHH8eabb7omkOQmJk+eLC1evLj3a4fDIUVEREjLly+/4PaLFy+Wrrvuuj73LVmyRJo2bdqQ5vREA9334eHh0iuvvNLnvttvv1266667hjSnpwMgffbZZxfd5rHHHpNSU1P73HfnnXdK8+bNG8Jk3uFy9v+3zZgxQ/rVr341ZHm8yUD3/TkpKSnSU0895fpAXmSw+/62226TfvzjH7skg1scUbHZbDh06BBmz57de59SqcTs2bOxZ8+eCz7m6quvxqFDh3pPURQWFuKrr77CggULhiWzpxjMvrdardDpdH3u8/X1xc6dO4c0KwF79uzp87MCgHnz5vX7syLyVE6nE21tbQgKChIdxescPnwYu3fvxowZM1zyfMJnpr0c9fX1cDgc581WGxYWhtOnT1/wMT/60Y9QX1+Pa665BpIkwW6344EHHuCpnwEazL6fN28eXnjhBUyfPh3x8fHYsmUL1q1bB4fDMRyRvVp1dfUFf1atra3o6uqCr6+voGREw+svf/kL2tvb8YMf/EB0FK8RFRWFuro62O12LFu2DPfff79LntctjqgMRmZmJp555hm89tpryM7Oxrp16/Dll1/ij3/8o+hoHu/FF1/EqFGjkJycDI1Gg4cffhj33nsvlEqP/XUjIhl5//338dRTT+Gjjz5CaGio6DheY8eOHTh48CBef/11rFixAh988IFLntctjqiEhIRApVKhpqamz/01NTWwWCwXfMz//u//4ic/+Ulvo0tLS0NHRwd+/vOf43/+53/4oXmZBrPvzWYz1q9fj+7ubjQ0NCAiIgK/+93vMHLkyOGI7NUsFssFf1YGg4FHU8grrF27Fvfffz8+/vjj806D0tCKi4sDcPbztqamBsuWLcMPf/jDK35et/i01mg0mDBhArZs2dJ7n9PpxJYtWzB16tQLPqazs/O8MqJSqQAAEpc3umyD2ffn6HQ6REZGwm6349NPP8Utt9wy1HG93tSpU/v8rABg06ZNl/xZEXmCDz74APfeey8++OAD3HjjjaLjeDWn0wmr1eqS53KLIyoAsGTJEixatAgTJ07E5MmTsWLFCnR0dODee+8FANx9992IjIzE8uXLAQALFy7ECy+8gIyMDEyZMgX5+fn43//9XyxcuLC3sNDlGei+37dvHyoqKjBu3DhUVFRg2bJlcDqdeOyxx0T+NdxSe3s78vPze78uKipCTk4OgoKCEBMTg6VLl6KiogJr1qwBADzwwAN45ZVX8Nhjj+G+++7D1q1b8dFHH+HLL78U9VdwawPd/wCQk5PT+9i6ujrk5ORAo9EgJSVluOO7tYHu+/fffx+LFi3Ciy++iClTpqC6uhrA2YH8RqNRyN/BXQ1037/66quIiYlBcnIygLOXiv/lL3/BL3/5S9cEcsm1Q8Pk5ZdflmJiYiSNRiNNnjxZ2rt3b+/3ZsyYIS1atKj3656eHmnZsmVSfHy8pNPppOjoaOmhhx6Smpqahj+4BxjIvs/MzJRGjx4tabVaKTg4WPrJT34iVVRUCEjt/rZt2yYBOO92bn8vWrRImjFjxnmPGTdunKTRaKSRI0dKq1atGvbcnmIw+/9C248YMWLYs7u7ge77GTNmXHR7unwD3fcvvfSSlJqaKvn5+UkGg0HKyMiQXnvtNcnhcLgkj0KSeB6EiIiI5MktxqgQERGRd2JRISIiItliUSEiIiLZYlEhIiIi2WJRISIiItliUSEiIiLZYlEhIiIi2WJRISIiItliUSEiIiLZYlEhIiIi2WJRISIiItliUSEiIiLZ+v+1WqWTEWNW0gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "df_pa = pd.DataFrame(ATE_estimates, columns=['Estimate'])\n",
    "g = sns.kdeplot(df_pa, fill=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outcome missing under nonignorable nonresponse\n",
    "Now consider a different setting, in which the outcomes are missing under nonignorable nonresponse assumptions in [Bia, Huber and Lafférs (2023)](https://doi.org/10.1080/07350015.2023.2271071). \n",
    "Let the covariance matrix $\\sigma^2_X$ again be such that $a_{ij} = 0.5^{|i - j|}$, but now $\\gamma_0 = 1$ and $\\sigma^2_{\\varepsilon, \\upsilon} = \\begin{pmatrix} 1 & 0.8 \\\\  0.8 & 1 \\end{pmatrix}$ to show a strong correlation between $\\varepsilon$ and $\\upsilon$. Let the vector of coefficients $\\beta$ again resemble a quadratic decay of coefficients importance; $\\beta_{0,j} = 0.4/j^2$ for $j = 1, \\ldots, p$.\n",
    "\n",
    "The directed acyclic graph (DAG) shows the the structure of the causal model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p align=\"center\">\n",
    "  <img src=\"figures/py_ssm.svg\" alt=\"Graph\">\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data\n",
    "\n",
    "We will again use the implemented data generating process `make_ssm_data` to generate data according to the simulation in Appendix E of [Bia, Huber and Lafférs (2023)](https://doi.org/10.1080/07350015.2023.2271071). We will again leave the default ATE equal to $\\theta_0=1$.\n",
    "\n",
    "In this setting, we are estimating the ATE under nonignorable nonresponse, so we set `mar=False`. Again, the selection indicator `S` can be set via `s_col`. Further, we need to specify an intrument via `z_col`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================== DoubleMLSSMData Object ==================\n",
      "Selection variable: s\n",
      "Outcome variable: y\n",
      "Treatment variable(s): ['d']\n",
      "Covariates: ['X1', 'X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8', 'X9', 'X10', 'X11', 'X12', 'X13', 'X14', 'X15', 'X16', 'X17', 'X18', 'X19', 'X20', 'X21', 'X22', 'X23', 'X24', 'X25', 'X26', 'X27', 'X28', 'X29', 'X30', 'X31', 'X32', 'X33', 'X34', 'X35', 'X36', 'X37', 'X38', 'X39', 'X40', 'X41', 'X42', 'X43', 'X44', 'X45', 'X46', 'X47', 'X48', 'X49', 'X50', 'X51', 'X52', 'X53', 'X54', 'X55', 'X56', 'X57', 'X58', 'X59', 'X60', 'X61', 'X62', 'X63', 'X64', 'X65', 'X66', 'X67', 'X68', 'X69', 'X70', 'X71', 'X72', 'X73', 'X74', 'X75', 'X76', 'X77', 'X78', 'X79', 'X80', 'X81', 'X82', 'X83', 'X84', 'X85', 'X86', 'X87', 'X88', 'X89', 'X90', 'X91', 'X92', 'X93', 'X94', 'X95', 'X96', 'X97', 'X98', 'X99', 'X100']\n",
      "Instrument variable(s): ['z']\n",
      "No. Observations: 8000\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "n_obs = 8000\n",
    "df = make_ssm_data(n_obs=n_obs, mar=False, return_type='DataFrame')\n",
    "dml_data = dml.DoubleMLSSMData(df, 'y', 'd', z_cols='z', s_col='s')\n",
    "print(dml_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimation\n",
    "\n",
    "We will again use the `DoubleMLSSM` class. \n",
    "\n",
    "Further,  will leave he learners for all nuisance functions to be the same as in the first setting, as the simulated quadratic decay of coefficients importance still holds.\n",
    "\n",
    "Now the learner `ml_g` is used to fit conditional expectations of the outcome $\\mathbb{E}[Y_i|D_i, S_i, X_i, \\Pi_i]$, whereas the learners `ml_m` and `ml_pi` will be used to estimate the treatment and selection propensity scores $P(D_i=1|X_i, \\Pi_i)$ and $P(S_i=1|D_i, X_i, Z_i)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml_g = LassoCV()\n",
    "ml_m = LogisticRegressionCV(penalty='l1', solver='liblinear')\n",
    "ml_pi = LogisticRegressionCV(penalty='l1', solver='liblinear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The score is now set to `'nonignorable'`, since the parameters of the DGP were set to satisfy the assumptions of outcomes missing under nonignorable nonresponse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================== DoubleMLSSM Object ==================\n",
      "\n",
      "------------------ Data Summary      ------------------\n",
      "Outcome variable: y\n",
      "Treatment variable(s): ['d']\n",
      "Covariates: ['X1', 'X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8', 'X9', 'X10', 'X11', 'X12', 'X13', 'X14', 'X15', 'X16', 'X17', 'X18', 'X19', 'X20', 'X21', 'X22', 'X23', 'X24', 'X25', 'X26', 'X27', 'X28', 'X29', 'X30', 'X31', 'X32', 'X33', 'X34', 'X35', 'X36', 'X37', 'X38', 'X39', 'X40', 'X41', 'X42', 'X43', 'X44', 'X45', 'X46', 'X47', 'X48', 'X49', 'X50', 'X51', 'X52', 'X53', 'X54', 'X55', 'X56', 'X57', 'X58', 'X59', 'X60', 'X61', 'X62', 'X63', 'X64', 'X65', 'X66', 'X67', 'X68', 'X69', 'X70', 'X71', 'X72', 'X73', 'X74', 'X75', 'X76', 'X77', 'X78', 'X79', 'X80', 'X81', 'X82', 'X83', 'X84', 'X85', 'X86', 'X87', 'X88', 'X89', 'X90', 'X91', 'X92', 'X93', 'X94', 'X95', 'X96', 'X97', 'X98', 'X99', 'X100']\n",
      "Instrument variable(s): ['z']\n",
      "No. Observations: 8000\n",
      "\n",
      "\n",
      "------------------ Score & Algorithm ------------------\n",
      "Score function: nonignorable\n",
      "\n",
      "------------------ Machine Learner   ------------------\n",
      "Learner ml_g: LassoCV()\n",
      "Learner ml_pi: LogisticRegressionCV(penalty='l1', solver='liblinear')\n",
      "Learner ml_m: LogisticRegressionCV(penalty='l1', solver='liblinear')\n",
      "Out-of-sample Performance:\n",
      "Regression:\n",
      "Learner ml_g_d0 RMSE: [[1.01990373]]\n",
      "Learner ml_g_d1 RMSE: [[1.19983954]]\n",
      "Classification:\n",
      "Learner ml_pi Log Loss: [[0.42934105]]\n",
      "Learner ml_m Log Loss: [[0.55828259]]\n",
      "\n",
      "------------------ Resampling        ------------------\n",
      "No. folds: 5\n",
      "No. repeated sample splits: 1\n",
      "\n",
      "------------------ Fit Summary       ------------------\n",
      "       coef   std err          t          P>|t|  2.5 %    97.5 %\n",
      "d  0.923517  0.037509  24.620995  7.528381e-134   0.85  0.997034\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "dml_ssm = dml.DoubleMLSSM(dml_data, ml_g, ml_m, ml_pi, score='nonignorable')\n",
    "dml_ssm.fit()\n",
    "\n",
    "print(dml_ssm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ATE estimates distribution\n",
    "\n",
    "Here we again add a small simulation where we generate multiple datasets, estimate the ATE and collect the results (this may take some time). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 0/200\n",
      "Iteration: 20/200\n",
      "Iteration: 40/200\n",
      "Iteration: 60/200\n",
      "Iteration: 80/200\n",
      "Iteration: 100/200\n",
      "Iteration: 120/200\n",
      "Iteration: 140/200\n",
      "Iteration: 160/200\n",
      "Iteration: 180/200\n"
     ]
    }
   ],
   "source": [
    "ml_g = LassoCV()\n",
    "ml_m = LogisticRegressionCV(penalty='l1', solver='liblinear')\n",
    "ml_pi = LogisticRegressionCV(penalty='l1', solver='liblinear')\n",
    "\n",
    "\n",
    "n_rep = 200\n",
    "ATE = 1.0\n",
    "\n",
    "ATE_estimates = np.full((n_rep), np.nan)\n",
    "\n",
    "np.random.seed(42)\n",
    "for i_rep in range(n_rep):\n",
    "    if (i_rep % int(n_rep/10)) == 0:\n",
    "        print(f'Iteration: {i_rep}/{n_rep}')\n",
    "    dml_data = make_ssm_data(n_obs=n_obs, mar=False)\n",
    "    dml_ssm = dml.DoubleMLSSM(dml_data, ml_g, ml_m, ml_pi, score='nonignorable')\n",
    "    dml_ssm.fit()\n",
    "\n",
    "    ATE_estimates[i_rep] = dml_ssm.coef.squeeze()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot the estimates distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAATNtJREFUeJzt3Xl8VOWhP/7P7JNktiwz2UhIyMK+uYCAiqh1QanVW+tVK4jXtha8lnKtlV5FbKtob6WoVaz9Wqi/itYFbeuCIgqoyE4A2bOQfU8mk8lMZjIz5/dHSEokLElm5jkz83m/XvNqM5zMfDhE5sNznuc5CkmSJBARERHJkFJ0ACIiIqIzYVEhIiIi2WJRISIiItliUSEiIiLZYlEhIiIi2WJRISIiItliUSEiIiLZYlEhIiIi2VKLDjAUgUAANTU1MBqNUCgUouMQERHReZAkCe3t7cjIyIBSefYxk4guKjU1NcjKyhIdg4iIiAahsrISw4YNO+sxEV1UjEYjgO7fqMlkEpyGiIiIzofD4UBWVlbv5/jZRHRR6bncYzKZWFSIiIgizPlM2+BkWiIiIpItoUXF7/fj0UcfRW5uLuLi4pCXl4ff/OY34A2diYiICBB86efpp5/GqlWr8Ne//hVjx47Frl27MH/+fJjNZjzwwAMioxEREZEMCC0qW7duxU033YQbbrgBAJCTk4PXX38dO3bsEBmLiIhkRpIk+Hw++P1+0VHoPKhUKqjV6qBsHSK0qEyfPh0vv/wyjh07hsLCQuzbtw9ffvklVqxY0e/xHo8HHo+n92uHwxGuqEREJIjX60VtbS1cLpfoKDQA8fHxSE9Ph1arHdLrCC0qDz/8MBwOB0aNGgWVSgW/348nnngCd955Z7/HL1++HI8//niYUxIRkSiBQABlZWVQqVTIyMiAVqvlBp8yJ0kSvF4vGhsbUVZWhoKCgnNu6nY2QovKm2++iddeew1r167F2LFjUVRUhEWLFiEjIwPz5s077fglS5Zg8eLFvV/3rMMmIqLo5PV6EQgEkJWVhfj4eNFx6DzFxcVBo9GgvLwcXq8Xer1+0K8ltKj84he/wMMPP4z//M//BACMHz8e5eXlWL58eb9FRafTQafThTsmEREJNpR/kZMYwfozE/on73K5TvuNqFQqBAIBQYmIiIhIToSOqMyZMwdPPPEEsrOzMXbsWOzduxcrVqzAPffcIzIWERHJXLXdjdYOb9jeLzFBi0xLXNjej/5NaFF5/vnn8eijj2LBggVoaGhARkYGfvKTn2Dp0qUiYxERkYxV29246plN6OwK3+i7XqPExv+5QmhZWbNmDRYtWgS73S4sgwhCi4rRaMTKlSuxcuVKkTGIiCiCtHZ40dkVwMJZ+WEpDtV2N174vBitHd4Bvd/dd9+Nv/71r6c9f+2112L9+vVn/d6cnBwsWrQIixYt6n3utttuw+zZs8/7/QdLboUoom9KSEREsSvTEofclATRMc7quuuuw+rVq/s8N9hFIXFxcYiLi73LT5xGTUREFCI6nQ5paWl9HomJiZAkCcuWLUN2djZ0Oh0yMjJ6bx1zxRVXoLy8HD//+c+hUCh6941Zs2YNLBZL72svW7YMkyZNwl/+8hdkZ2fDYDBgwYIF8Pv9+N3vfoe0tDTYbDY88cQTfTKtWLEC48ePR0JCArKysrBgwQI4nU4AwKZNmzB//ny0tbX1vveyZcsAdG+6+uCDDyIzMxMJCQmYOnUqNm3aFPJzyBEVilqBgISPvqnDpqMN2Fthx3cnZeC/r8znZlFEJNw777yDP/zhD3jjjTcwduxY1NXVYd++fQCAdevWYeLEifjxj3+MH/3oR2d9nZKSEnz00UdYv349SkpK8P3vfx+lpaUoLCzE5s2bsXXrVtxzzz24+uqrMXXqVADdy4afe+455ObmorS0FAsWLMBDDz2EF198EdOnT8fKlSuxdOlSHD16FABgMBgAAPfffz8OHTqEN954AxkZGXj33Xdx3XXX4cCBAygoKAjZuWJRoaj1/74sxZMfHkF2UjwyLHqs2HAMxQ3t+N33J0KvUYmOR0Qx4P333+/9oO/xq1/9Cnq9Hmlpabj66quh0WiQnZ2NKVOmAACSkpKgUqlgNBqRlpZ21tcPBAL4y1/+AqPRiDFjxmDWrFk4evQoPvzwQyiVSowcORJPP/00Pv/8896icuq8l5ycHPz2t7/FfffdhxdffBFarRZmsxkKhaLPe1dUVGD16tWoqKhARkYGAODBBx/E+vXrsXr1ajz55JPBOF39YlGhqHSoxoHfrT+KGyek486pwwEAl+Q248VNJfD49uJPd10kOCERxYJZs2Zh1apVfZ5LSkpCR0cHVq5ciREjRuC6667D7NmzMWfOHKjVA/tYzsnJgdFo7P06NTUVKpWqzx5lqampaGho6P36008/xfLly3HkyBE4HA74fD50dnbC5XKdcfffAwcOwO/3o7CwsM/zHo8HycnJA8o8UCwqFHU6u/z42Rt7MSwxDj+46N+3WJg6Ihm+gIQ/fl6MnSdacHFOksCURBQLEhISkJ+ff9rzSUlJOHr0KD799FNs2LABCxYswP/93/9h8+bN0Gg05/363z5WoVD0+1zPRqonTpzAjTfeiJ/+9Kd44oknkJSUhC+//BL/9V//Ba/Xe8ai4nQ6oVKpsHv3bqhUfUekvz1iFGycTEtR55Uvy3CiuQMLrsiHRtX3R3xaXjJyUxLw1EdHIEmSoIRERN2reObMmYPnnnsOmzZtwtdff40DBw4AALRaLfx+f9Dfc/fu3QgEAnjmmWdwySWXoLCwEDU1NX2O6e+9J0+eDL/fj4aGBuTn5/d5nOvy1FBxRIWiSiAgYe32CkzPS0FW0un/MlAqFLjtoiw8tf4IPjvSgKtGpwpISUTBUG13y/59PB4P6urq+jynVqvx/vvvw+/3Y+rUqYiPj8ff/vY3xMXFYfjw7kvVOTk52LJlC/7zP/8TOp0OKSkpQ/o99MjPz0dXVxeef/55zJkzB1999RVeeumlPsfk5OTA6XRi48aNmDhxIuLj41FYWIg777wTc+fOxTPPPIPJkyejsbERGzduxIQJE3DDDTcEJV9/WFQoqmwtaUa13Y0fXTbijMdMGGbG2AwTfrf+KK4cZeMqIKIIk5ighV6jxAufF4ftPfUaJRITtAP+vvXr1yM9Pb3PcyNHjsRTTz2Fp556CosXL4bf78f48ePxr3/9q3e+x69//Wv85Cc/QV5eHjweT9BGgCdOnIgVK1bg6aefxpIlS3D55Zdj+fLlmDt3bu8x06dPx3333YfbbrsNzc3NeOyxx7Bs2TKsXr0av/3tb/E///M/qK6uRkpKCi655BLceOONQcl2Jgopgse/HQ4HzGYz2traYDKZRMchGfjvtXuwu8KO339/wlkLyDfVbXjiw8NYt2A6LshODGNCIhqIzs5OlJWVITc3F3q9vvd53utH/s70ZwcM7PObIyoUNewuLz4+WI9bLxp2zlGSMekmJCdo8d7eahYVogiUaYljcYgRnExLUeO9vdXwSxIuK7Ce81ilUoFpecn4574adPnDd2MzIiIaGBYVihrvFdVgcpYF5rjzW9p3aX4K7K4ubDnWGOJkREQ0WCwqFBXaXF3YX2XHpGzLeX/P8OQEZCfF49291aELRkREQ8KiQlHh69ImBCRgQqZ5QN83Iz8FGw7Vo72zK0TJiCgYInjdR8wK1p8ZiwpFhS+Lm5Bu1sNq1J/74FPMyEuGxxfA50d5+YdIjnp2WXW5XIKT0ED1/JkNZKfd/nDVD0WFL441YWzGwEZTACDZoEN2Ujw2H23EdydmhCAZEQ2FSqWCxWLpvVdNfHw89z6SOUmS4HK50NDQAIvFctqW+wPFokIRr7LFhfIWF/7jgmGD+v4Jw8zYdKwBgYAEpZJ/ARLJTc8W7afeWI/kz2KxBGV7fRYVinhfFjdBqQDGZAxu079JWRa8v78Wh2odGDfAOS5EFHoKhQLp6emw2Wzo6uJ8skig0WiGPJLSg0WFIt4XxxuRZzMgQTe4H+eRqUboNUpsPtbIokIkYyqVKmgffhQ5OJmWIlogIOGr4maMH8T8lB5qlRJjM8zYdJTDykREcsOiQhGtrLkDbe4ujEwzDul1Jg6zYE+5HQ4uUyYikhUWFYpo+yrtAIARVsOQXmfiMDP8koStxU1BSEVERMHCokIRbX9VGzIsehgGOT+lh82kR4ZFjy+Os6gQEckJiwpFtL2VrchNGdpoSo9RaSZsL2sJymsREVFwsKhQxPL6Ajhc0458a0JQXm9UmhHFDU40Oz1BeT0iIho6FhWKWEfr2uH1B5A3xPkpPcakd+/DsvNEa1Bej4iIho5FhSJWUZUdKqUCw5ODM6KSbNDBZtRhBy//EBHJBosKRaz9lXYMT4qHVh28H+ORaUZsK20O2usREdHQsKhQxCqqtCM3JTijKT1Gp5lwuNbB/VSIiGSCRYUiktPjQ3GDE3m24MxP6TEq3QgJwG7OUyEikgUWFYpIB6vbIAFBm0jbI82kR1K8hsuUiYhkgkWFItKRunaoVQpkWPRBfV2FQoGRaSZsL+M8FSIiOWBRoYh0pK4dwyxxUCuD/yNcmGrAN9Vt8Pj8QX9tIiIaGKFFJScnBwqF4rTHwoULRcaiCHCkzoHMxPiQvHZBqhFdfgkHaxwheX0iIjp/QovKzp07UVtb2/vYsGEDAODWW28VGYtkTpIkHKtrR3ZiXEhef3hSPLQqJfZW2EPy+kREdP6EFhWr1Yq0tLTex/vvv4+8vDzMnDlTZCySuWq7Gx1eP4YlhWZERa1SIteagD0VXPlDRCTa0G45G0Rerxd/+9vfsHjxYigUin6P8Xg88Hj+fR8Wh4ND87HoWH07ACA7REUFAApsBuwuZ1EhIhJNNpNp33vvPdjtdtx9991nPGb58uUwm829j6ysrPAFJNk4UteOeK0KyQnakL1Hvs2A2rZO1Ds6Q/YeRER0brIpKq+88gquv/56ZGRknPGYJUuWoK2trfdRWVkZxoQkF8fq2pGVGH/GkbdgKLAZAQB7OKpCRCSULC79lJeX49NPP8W6devOepxOp4NOpwtTKpKrw3XtGBaiibQ9khK0sBp02Ftpx/Xj00P6XkREdGayGFFZvXo1bDYbbrjhBtFRSOa6/AGUNDiRFcL5KT3yOU+FiEg44UUlEAhg9erVmDdvHtRqWQzwkIydaOqALyCFragcqGqD1xcI+XsREVH/hBeVTz/9FBUVFbjnnntER6EIcKSue8VPVogv/QDdK3+8/gCOnnxPIiIKP+FDGNdccw0kSRIdgyLEsfp2JMZrYNRrQv5ew5MToFIqUFRlx/hh5pC/HxERnU74iArRQJQ0OpFhCf1oCgBo1UrkJMdjX6U9LO9HRESnY1GhiFLc4ES6Obh3TD6b3BQDiriVPhGRMCwqFDECAQnlzS6km8MzogIA+bYElDQ60d7ZFbb3JCKif2NRoYhRbXfD4wsgwxK+EZU8qwESgAPVbWF7TyIi+jcWFYoYZU0dABDWEZUMcxziNCrsq2RRISISgUWFIkZpoxMalQJWQ/h2J1YqFRhhTcC+Sm78RkQkAosKRYzSpg6kmfVQKkN3j5/+5FkNKOKIChGRECwqFDFKGp1IM4VvfkqPPKsBdQ7eSZmISAQWFYoYJY0dYZ2f0iPPmgAAKOJ+KkREYceiQhHB5fWhrq0zrCt+eiQlaJEYr8H+KnvY35uIKNaxqFBEELHip4dCocAIq4EjKkREArCoUEToKSoZAooKAIxIScD+qjYEArwvFRFROLGoUEQobeyASa+GQS/mPpr5NgPaO3040dwh5P2JiGIViwpFhNJGJ9LDdDPC/oxIMQAA9ldxmTIRUTixqFBEKGnsELI0uYdBr0a6Wc95KkREYcaiQhGhvLl7szeRRqQkYB9X/hARhRWLCsme3eWFo9OHVKPgomI14FCNA13+gNAcRESxhEWFZK+82QUASDWF7x4//cm3GeDxBXC0rl1oDiKiWMKiQrJX3tJTVMSOqOQkJ0CpAC//EBGFEYsKyV5FcweMejUSdGKWJvfQqpUYnhyPfZxQS0QUNiwqJHvlzS7hoyk9RqQYsLfCLjoGEVHMYFEh2TvR3AGbUez8lB55NgOKG5xwenyioxARxQQWFZI9OY2o5FsNkAAc4MZvRERhwaJCstbZ5UdDu0f4ip8emZY4xGlU3PiNiChMWFRI1ip6VvwI3kOlh1KpwAhrAooqW0VHISKKCSwqJGsnTt41OVXwrrSnyrMaOKJCRBQmLCokaxUtLujUSljiNKKj9Mq3GlDv8KCurVN0FCKiqMeiQrLWM5FWoVCIjtIrz9Z9J2Ve/iEiCj0WFZK1chktTe6RlKBFcoIWRZVc+UNEFGosKiRrJ2S0NPlUeVYD9lZwRIWIKNRYVEi2fP4Aauxu2SxNPlWezYD9VW3w8U7KREQhxaJCslXb1glfQIJNJkuTT1WYaoC7y48jvJMyEVFIsaiQbFWe3ENFbnNUgO57/qiUCl7+ISIKMRYVkq2qVjcUAFJkWFS0aiVyk+OxhzcoJCIKKeFFpbq6Gj/84Q+RnJyMuLg4jB8/Hrt27RIdi2SgstWFpAQtNCrhP6b9yrcZsbucIypERKEk9BOgtbUVM2bMgEajwUcffYRDhw7hmWeeQWJioshYJBOVLS5Zjqb0KEw1oKLFhSanR3QUIqKopRb55k8//TSysrKwevXq3udyc3MFJiI5qWhxwWqQb1EpSDUCAPaUt+KasWmC0xARRSehIyr//Oc/cdFFF+HWW2+FzWbD5MmT8ec//1lkJJKRyla3LCfS9khO0CIpQct5KkREISS0qJSWlmLVqlUoKCjAxx9/jJ/+9Kd44IEH8Ne//rXf4z0eDxwOR58HRafOLj8a2z2wyrioKBQK5NsM2F3eIjoKEVHUElpUAoEALrjgAjz55JOYPHkyfvzjH+NHP/oRXnrppX6PX758Ocxmc+8jKysrzIkpXKrtbgDyXJp8qgKbAQeq2tDFjd+IiEJCaFFJT0/HmDFj+jw3evRoVFRU9Hv8kiVL0NbW1vuorKwMR0wSoGcPFTmPqADAyFQjOn0BfFPN+/4QEYWC0Mm0M2bMwNGjR/s8d+zYMQwfPrzf43U6HXQ6eX9wUXBUtrqhVABJCfL+885NSYBOrcTOEy2YnM3VakREwSZ0ROXnP/85tm3bhieffBLFxcVYu3YtXn75ZSxcuFBkLJKBqlYXrEYdVEqF6ChnpVYpkW8zYEcZ56kQEYWC0KJy8cUX491338Xrr7+OcePG4Te/+Q1WrlyJO++8U2QskoGqFreslyafamSaETtPtCIQkERHISKKOkIv/QDAjTfeiBtvvFF0DJKZ8haX7Oen9BidZsK6PdUobnSi8OTeKkREFBzy3JucYl5ViwspETKikm/rvkHhdl7+ISIKOhYVkh2nxwe7uws2k150lPOi16iQm5KAnSwqRERBx6JCslPV2r00We57qJxqVJoRO8paIEmcp0JEFEwsKiQ7lS3dm71FyqUfoHtCbZ2jE1WtbtFRiIiiCosKyU5VqwsalQKWeI3oKOdtVKoJADhPhYgoyFhUSHaqW91IMeigVMh7D5VTGfRq5CTH4+uSZtFRiIiiCosKyU613R1Rl316jMkw46viJs5TISIKIhYVkp3KVjdSDFrRMQZsXIYJdY5OlDV1iI5CRBQ1WFRIdqpbI2cPlVONSjNBpVRgKy//EBEFDYsKyYrL60OrqytidqU9VZxWhXyrAV8VN4mOQkQUNVhUSFZq7JG3NPlUYzJM+Lqkmff9ISIKEhYVkpXKk/uQROKICtA9T8Xu7sLhOofoKEREUYFFhWSlutUNpQJIjI+8ybQAUJBqhFalxNZizlMhIgoGFhWSlZ6lySpl5OyhciqNSolRaUZ8cbxRdBQioqjAokKy0rPZWySbMMyCbaUtcHv9oqMQEUU8FhWSlcpWF5IjcA+VU03KssDrD2BbKS//EBENFYsKyUp1qxvWCB9RybDoYTPqsOlog+goREQRj0WFZMPj86Ox3RPxl34UCgUmDLPg86Ocp0JENFQsKiQbtfZOSIjcpcmnmpRlQUWLi9vpExENEYsKyUZ1hG/2dqqxGSZoVApe/iEiGiIWFZKN6pObvUX6ZFoA0GtUGJVmwudHWFSIiIaCRYVko8ruRlKCFhpVdPxYTsrqXqbc4fGJjkJEFLGi4xOBokL3HiqRP5rS46LhifD6A9h8jJNqiYgGi0WFZKOq1YXkhMifn9LDZtIjJzkeHx+sEx2FiChisaiQbHRvnx89IyoAcOHwRHx2uAFd/oDoKEREEYlFhWTBH5BQ19YZFSt+TnVRThLaPT5sL20RHYWIKCKxqJAsNLZ74AtIUVdUhifFw2rU8fIPEdEgsaiQLPTuoRIFm72dSqFQ4MLhifjkUB0CAUl0HCKiiMOiQrLw783eomuOCgBcPDwR9Q4PiqrsoqMQEUUcFhWShRq7GwlaFeK1atFRgm5UmgmJ8Rp8sL9WdBQioojDokKyUN3qRnKUzU/poVQqcHFOEt7fX8PLP0REA8SiQrIQjUuTTzUtLxn1Dg92lbeKjkJEFFFYVEgWqlpdUTuiAgCFqUYkJ2jx/v4a0VGIiCIKiwrJQq09+vZQOZVSocDU3CR8cKAWfl7+ISI6b0KLyrJly6BQKPo8Ro0aJTISCeDo7EK7xxfVl36A7ss/zU4vtpU2i45CRBQxhC+xGDt2LD799NPer9Vq4ZEozGp6lyZH74gKAORZDUg16fCPomrMyE8RHYeIKCIIv/SjVquRlpbW+0hJ4V/gsaa6NTaKikKhwIy8FHx4oA6dXX7RcYiIIoLwonL8+HFkZGRgxIgRuPPOO1FRUSE6EoVZjd0NtVIBS7xGdJSQu7QgBU6PD58erhcdhYgoIggtKlOnTsWaNWuwfv16rFq1CmVlZbjsssvQ3t7e7/EejwcOh6PPgyJfld2NZIMWSoVCdJSQSzfHocBmwLo91aKjEBFFBKFF5frrr8ett96KCRMm4Nprr8WHH34Iu92ON998s9/jly9fDrPZ3PvIysoKc2IKhZooX/HzbTPyU7D5WCOanR7RUYiIZE/4pZ9TWSwWFBYWori4uN9fX7JkCdra2noflZWVYU5IoVDV6kJSQnSv+DnVtLxkAMC/9nFPFSKic5FVUXE6nSgpKUF6enq/v67T6WAymfo8KPLV2N0xNaJi0mswaZgF7/DyDxHROQktKg8++CA2b96MEydOYOvWrbj55puhUqlw++23i4xFYeT1BdDg8MRUUQGAywutOFDdhmP1/c/HIiKibkKLSlVVFW6//XaMHDkSP/jBD5CcnIxt27bBarWKjEVhVO/ohARE/WZv33ZBtgVGvRrv7K4SHYWISNaE7q72xhtviHx7koHqGNns7dvUKiWm56XgnT1V+MW1I6FWyeoqLBGRbPBvRxKqZ7O35BgbUQGAmYVWNDm9+OJ4k+goRESyxaJCQtXY3TDHaaBTq0RHCbuc5HhkJ8Xjrd1cvUZEdCYsKiRUTZs75uan9FAoFLi8wIoNh+phd3lFxyEikiUWFRKqqtUdU3uofNuM/GT4AxL+tb9WdBQiIlliUSGhqltjaw+Vb7PEazEpy4K3dvHyDxFRf1hUSBhJkk5e+ondogJ076myv6oNx7mnChHRaVhUSJhWVxc6uwIxueLnVBdkJ8KgU+PtPdxThYjo21hUSJiaGN1D5ds0KiWm5yXj3T3V8Ack0XGIiGSFRYWEqWplUelxeaEVDe0efHG8UXQUIiJZYVEhYWrsbmhVSpj0QjdIloURKQnItMTh3b28USER0alYVEiYGrsbKUYtFAqF6CjCKRQKXJqfgo8P1qHD4xMdh4hINgZVVEpLS4Odg2JQtd2N5ARe9ukxIz8FnV0BrP+mTnQUIiLZGFRRyc/Px6xZs/C3v/0NnZ2dwc5EMaKqNXZ3pe2P1ajD6HQjL/8QEZ1iUEVlz549mDBhAhYvXoy0tDT85Cc/wY4dO4KdjaJcTZsbyZxI28el+VZsLWlCvYP/ACAiAgZZVCZNmoRnn30WNTU1+Mtf/oLa2lpceumlGDduHFasWIHGRq5coLPr7PKj2enliMq3TM1NgkqpwD+LakRHISKShSFNplWr1bjlllvw1ltv4emnn0ZxcTEefPBBZGVlYe7cuait5f1LqH+1bd0jBpyj0leCTo3JWYn45z4WFSIiYIhFZdeuXViwYAHS09OxYsUKPPjggygpKcGGDRtQU1ODm266KVg5KcpUn9xDxWpkUfm2S0Yk40B1G8qbO0RHISISblAbWKxYsQKrV6/G0aNHMXv2bLz66quYPXs2lMru3pObm4s1a9YgJycnmFkpilTbXVAASI7hOyefyeRsC/QaJd7fX4uFs/JFxyEiEmpQIyqrVq3CHXfcgfLycrz33nu48cYbe0tKD5vNhldeeSUoISn6VNs7kZighVrFrXy+Ta9R4YJsXv4hIgIGOaKyYcMGZGdnn1ZOJElCZWUlsrOzodVqMW/evKCEpOhTzaXJZ3XJiGSs2HAMxQ3tyLcZRcchIhJmUP+czcvLQ1NT02nPt7S0IDc3d8ihKPpVtbo4kfYsJg6zIF6rwr/2cUI6EcW2QRUVSer/Dq9OpxN6vX5IgSg2VNs5onI2WrUSF2Yn4oMDLCpEFNsGdOln8eLFALrvS7J06VLEx8f3/prf78f27dsxadKkoAak6OMPSKhr60TKGI6onM3FuUlYseEYShqdyLMaRMchIhJiQEVl7969ALpHVA4cOACt9t//ItZqtZg4cSIefPDB4CakqNPY7oEvICGFl37OasIwM3RqJT4+WIcFV3D1DxHFpgEVlc8//xwAMH/+fDz77LMwmUwhCUXRrdruAgCkcA+Vs9KpVZg4zIL137CoEFHsGtQcldWrV7Ok0KBVndzsjXNUzu2inETsr2pDXRvv/UNEsem8R1RuueUWrFmzBiaTCbfccstZj123bt2Qg1H0qrF3IkGnQrx2UKvjY8rk7ESolAp8cqgOc6fliI5DRBR25/1JYTaboVAoev8/0WBV211I4V2Tz4tBp8bYdBPWf8OiQkSx6byLyurVq/v9/0QDVd3q5kTaAbgoJxF/3VqONlcXzPEa0XGIiMJqUHNU3G43XC5X79fl5eVYuXIlPvnkk6AFo+hV1epGMuennLcLshPhlyRsOd4oOgoRUdgNqqjcdNNNePXVVwEAdrsdU6ZMwTPPPIObbroJq1atCmpAii6SJJ3c7I0jKucr2aDD8OR4fHakQXQUIqKwG1RR2bNnDy677DIAwNtvv420tDSUl5fj1VdfxXPPPRfUgBRdHG4fXF4/i8oATcqy4POjDfAH+t8VmogoWg2qqLhcLhiN3TdK++STT3DLLbdAqVTikksuQXl5eVADUnSpOrmHitXISz8DMTkrEXZXF4oq7aKjEBGF1aCKSn5+Pt577z1UVlbi448/xjXXXAMAaGho4P4qdFY19u79QDiiMjAFNgOMejU+5+UfIooxgyoqS5cuxYMPPoicnBxMnToV06ZNA9A9ujJ58uSgBqToUt3qgkalgCmOq1cGQqlUYEKmGRuP1IuOQkQUVoMqKt///vdRUVGBXbt2Yf369b3PX3XVVfjDH/4wqCBPPfUUFAoFFi1aNKjvp8hQ1do9kVZ5ck8eOn+TsxNxuLadu9QSUUwZ9NagaWlpSEtL6/PclClTBvVaO3fuxJ/+9CdMmDBhsHEoQlS1umHlZZ9BmTDMDAWAzccacNvF2aLjEBGFxaBGVDo6OvDoo49i+vTpyM/Px4gRI/o8BsLpdOLOO+/En//8ZyQmJg4mDkWQqlYXb0Y4SEa9BnnWBHxxvEl0FCKisBnUiMq9996LzZs346677kJ6enrv1vqDsXDhQtxwww24+uqr8dvf/vasx3o8Hng8nt6vHQ7HoN+XxKhqdWNMBm/BMFjjMs3YfKwRgYAEpZKXz4go+g2qqHz00Uf44IMPMGPGjCG9+RtvvIE9e/Zg586d53X88uXL8fjjjw/pPUmcDo8PdncXrBxRGbTxwyx4r6gGh2odGJfJwkdE0W9Ql34SExORlJQ0pDeurKzEz372M7z22mvQ6/Xn9T1LlixBW1tb76OysnJIGSi8qu1uAOAclSEotBmg1yh5+YeIYsagispvfvMbLF26tM/9fgZq9+7daGhowAUXXAC1Wg21Wo3Nmzfjueeeg1qtht/vP+17dDodTCZTnwdFjqrWns3eWFQGS61SYnSaCVuO8b4/RBQbBnXp55lnnkFJSQlSU1ORk5MDjabvnhh79uw552tcddVVOHDgQJ/n5s+fj1GjRuGXv/wlVCrVYKKRjFW1uqFWKmDhHYCHZPwwM17fUQG31484Lf87IaLoNqii8r3vfW/Ib2w0GjFu3Lg+zyUkJCA5Ofm05yk6VHMPlaAYn2nGq34JO060YGahVXQcIqKQGlRReeyxx4Kdg2JAVasbKbzHz5BlWuKQnKDFl8cbWVSIKOoNesM3u92Ot99+GyUlJfjFL36BpKQk7NmzB6mpqcjMzBzUa27atGmwcSgCVLS4kJLA+SlDpVAoMDrdhK0lzaKjEBGF3KAm0+7fvx+FhYV4+umn8fvf/x52ux0AsG7dOixZsiSY+SiKVNvdnEgbJGMyTDhU40Cbq0t0FCKikBpUUVm8eDHuvvtuHD9+vM/S4tmzZ2PLli1BC0fRw+31o6XDy6ISJGPTTZAAbC/jqAoRRbdBFZWdO3fiJz/5yWnPZ2Zmoq6ubsihKPpU27uXJqdwD5WgsJn0sBp12FbaIjoKEVFIDaqo6HS6frevP3bsGKxWTu6j01W1ntzsjSMqQTMm3YStJdz4jYii26CKyne/+138+te/RldX9/VxhUKBiooK/PKXv8R//Md/BDUgRYeqVjdUSgUS47nqJ1hGp5twpK4drR1e0VGIiEJmUEXlmWeegdPphNVqhdvtxsyZM5Gfnw+j0Ygnnngi2BkpClS1upFi0ELFG+kFzdiM7p2Zt5fx8g8RRa9BLU82m83YsGEDvvrqK+zbtw9OpxMXXHABrr766mDnoyhRbXdzfkqQpRh0SDXpsK20GdeNSxMdh4goJAZcVAKBANasWYN169bhxIkTUCgUyM3NRVpaGiRJgoK7jlI/yps7WFRCgPNUiCjaDejSjyRJ+O53v4t7770X1dXVGD9+PMaOHYvy8nLcfffduPnmm0OVkyJcVYsLNk6kDbpRaSYcq3fC7uI8FSKKTgMaUVmzZg22bNmCjRs3YtasWX1+7bPPPsP3vvc9vPrqq5g7d25QQ1Jk6/D40OLq4oqfEBiVZgQA7DrRiqvHpApOQ0QUfAMaUXn99dfxq1/96rSSAgBXXnklHn74Ybz22mtBC0fRoWdpss2oP8eRNFBWow7JCVrsPMEJtUQUnQZUVPbv34/rrrvujL9+/fXXY9++fUMORdGlsqV7szeOqASfQqFAYZoR27hDLRFFqQEVlZaWFqSmnnl4OTU1Fa2trUMORdGlstUFjUoBS7xGdJSoNDrNiIPVDri8PtFRiIiCbkBFxe/3Q60+87QWlUoFn49/WVJflS1u2Ix6KLkiLCRGpZngC0goqrCLjkJEFHQDmkwrSRLuvvtu6HT9D+F7PJ6ghKLoUtHSwcs+IZSZGAeDTo0dJ1owPT9FdBwioqAaUFGZN2/eOY/hih/6tooWF4YnJ4iOEbWUCgVGphqxgzvUElEUGlBRWb16dahyUJSSJAlVrW5cnJMkOkpUG5lmxLq9VejyB6BRDerOGEREssS/0SikWl1dcHn9vPQTYqPSjOjsCuBQzel3NSciimQsKhRSPUuTuYdKaOWmJECjUmBXOVfdEVF0YVGhkKps5R4q4aBWKZFnNWAXN34joijDokIhVdniRoJOBYNuUDfqpgEoTDVid3krJEkSHYWIKGhYVCikKlpcvOwTJoWpRjS0e3pvWUBEFA1YVCikKltcsBp42SccClINAIA9FZynQkTRg0WFQqqixcX5KWFi0muQaYnDrhMsKkQUPVhUKGT8AQk1djdsJhaVcMm3GXgnZSKKKiwqFDI1djd8AQmpnKMSNiNTjThW3472zi7RUYiIgoJFhUKm4uQeKqkmFpVwKUw1IiABRZV20VGIiIKCRYVC5kRzB5QKIMWoFR0lZmRY9DDo1NjNjd+IKEqwqFDIVDS7YDXqoVbyxyxcFAoFCmwG7GFRIaIowU8QCpnyZhdSueIn7PJtBuytsCMQ4MZvRBT5WFQoZMqaO7jiR4DCVCPaPT6UNDpFRyEiGjIWFQoJSZJQ0eziRFoB8qwGKBXc+I2IogOLCoVEk9MLd5efRUWAOK0KWUnxnFBLRFGBRYVCoqKlAwCXJotSYDOwqBBRVBBaVFatWoUJEybAZDLBZDJh2rRp+Oijj0RGoiA50dS9h4qNk2mFKLAZUdLYgTYXN34josgmtKgMGzYMTz31FHbv3o1du3bhyiuvxE033YSDBw+KjEVBUN7iQmK8BnqNSnSUmNRzg8K9lRxVIaLIJrSozJkzB7Nnz0ZBQQEKCwvxxBNPwGAwYNu2bSJjURBUNHfwso9AaSY9THo19lTYRUchIhoStegAPfx+P9566y10dHRg2rRp/R7j8Xjg8Xh6v3Y4HOGKRwN0otnFyz4CKRQK5HPjNyKKAsIn0x44cAAGgwE6nQ733Xcf3n33XYwZM6bfY5cvXw6z2dz7yMrKCnNaOl/lHFERrsBmRFElN34josgmvKiMHDkSRUVF2L59O376059i3rx5OHToUL/HLlmyBG1tbb2PysrKMKel89He2YVWVxeLimAFqQY4PT4cb+DGb0QUuYRf+tFqtcjPzwcAXHjhhdi5cyeeffZZ/OlPfzrtWJ1OB52OlxPkrry5567J/LMSqWfjt70VrRiZZhQdh4hoUISPqHxbIBDoMw+FIk9ZU/ceKmmmOMFJYpteo0J2Ujx3qCWiiCZ0RGXJkiW4/vrrkZ2djfb2dqxduxabNm3Cxx9/LDIWDVFZUwdMejUMeuEDdjEvnxu/EVGEE/pJ0tDQgLlz56K2thZmsxkTJkzAxx9/jO985zsiY9EQlTV1IN3C0RQ5KLAZ8enhBrS5u2CO04iOQ0Q0YEKLyiuvvCLy7SlEShqdSONEWlkosHVv/FZUacfMQqvgNEREAye7OSoU2SRJ6h5RMbOoyEGaWQ+jXs39VIgoYrGoUFA1d3jR3ulDupmXfuSgd+M3TqglogjFokJB1bPihyMq8pFvNaCoghu/EVFkYlGhoCpr7IAC4GZvMlKYakS7x4fiRm78RkSRh0WFgqq0qQNWow5aNX+05CLf1r3xG+epEFEk4qcJBVVZE1f8yA03fiOiSMaiQkFV2tiBNM5PkR1u/EZEkYpFhYLGH5BQ3uziRFoZKrAZUdLYAbvLKzoKEdGAsKhQ0NTY3fD6A0jj0mTZKUztvinh3kq72CBERAPEokJBw6XJ8pVq0sGkV2MvL/8QUYRhUaGgKW10Qq1UwGrQiY5C36JQKFBgM2I3J9QSUYRhUaGgKW50IsMSB6VSIToK9SM/tXvjNz83fiOiCMKiQkFzvN7Jyz4yVphqRIfXj2P17aKjEBGdNxYVCprjDU4MS+REWrnKsyZApVRwmTIRRRQWFQqKlg4vWjq8yLTEi45CZ6BTq5CbHM+iQkQRhUWFgqK4ofs+MhxRkbf8VCN2nmgRHYOI6LyxqFBQFDc4oVSAu9LK3MhUI6pa3WhwdIqOQkR0XlhUKCiON7QjzayHRsUfKTnr2fiNl3+IKFLwU4WC4ni9E5kWXvaRu6QELaxGHXaxqBBRhGBRoaA43tDOohIhCm0G7CrnPBUiigwsKjRk7Z1dqHd4kJnIFT+RoDDNiIPVDnR2+UVHISI6JxYVGrKeFT8cUYkMhalG+AIS9vEGhUQUAVhUaMiONzihAJBh4YqfSJCdGI84jYrzVIgoIrCo0JCVNDhhM+mgU6tER6HzoFQqUJhqwI4yzlMhIvljUaEhO1bfjgwzL/tEkpFpJuwub+UNColI9lhUaMgO17YjK4kTaSPJqDQjnB4fjtQ5REchIjorFhUaErvLizpHJ4Yns6hEkjyrARqVAjt5+YeIZI5FhYbkcG07ACCbIyoRRatWYoTVgB287w8RyRyLCg3JkToHNCoF0jlHJeKMTDViR1kLJInzVIhIvlhUaEgO1zqQlRgPlVIhOgoN0Mg0I5qcXpQ3u0RHISI6IxYVGpJDNQ5OpI1QI1ONUAC8/ENEssaiQoPm8wdwvMHJ+SkRKkGnRnZyPCfUEpGssajQoJ1o7oDHF+CKnwg2Ks2EbaXNomMQEZ0RiwoNGlf8RL4x6SZUtrpRbXeLjkJE1C+hRWX58uW4+OKLYTQaYbPZ8L3vfQ9Hjx4VGYkG4HCtA8kJWhj1GtFRaJBGpxsBANtKOKpCRPIktKhs3rwZCxcuxLZt27BhwwZ0dXXhmmuuQUdHh8hYdJ6O1Do4mhLhjHoNhifH8/IPEcmWWuSbr1+/vs/Xa9asgc1mw+7du3H55ZcLSkXn61BtO6bkJomOQUM0Os2Er1lUiEimZDVHpa2tDQCQlNT/h5/H44HD4ejzIDGanR7UOTqRw4m0EW9MhglVrW5UtXI/FSKSH9kUlUAggEWLFmHGjBkYN25cv8csX74cZrO595GVlRXmlNRjf3V3qRxhNQhOQkM1Os0EBYBtpVymTETyI5uisnDhQnzzzTd44403znjMkiVL0NbW1vuorKwMY0I61YGqNhj1atiMOtFRaIgMejXnqRCRbAmdo9Lj/vvvx/vvv48tW7Zg2LBhZzxOp9NBp+MHoxzsq7IjNyUBCgW3zo8Go9JN2FrSJDoGEdFphI6oSJKE+++/H++++y4+++wz5ObmioxD50mSJOyrtGNESoLoKBQkYzNMqLF3oryZK+6ISF6EjqgsXLgQa9euxT/+8Q8YjUbU1dUBAMxmM+LieDdeuap3eNDk9HJ+ShQZk26CUgF8WdyE4cksoEQkH0JHVFatWoW2tjZcccUVSE9P7338/e9/FxmLzmFflR0AOKISReK1ahTYjPjiOC//EJG8CB1RkSRJ5NvTIB2oakNivAZJCVrRUSiIxmWasOFQPfwBCSol5x4RkTzIZtUPRQ5OpI1O4zLNcHT68M3JpedERHLAokIDIkkS9le1cX5KFMq3GRCnUeHLYl7+ISL5YFGhAalscaPN3YU8K+enRBu1Uokx6SbOUyEiWWFRoQHZW9kKABiRwhGVaDQu04Td5S1we/2ioxARAWBRoQHaeaIFmZY4mOI0oqNQCIzPtKDLL2FbGXepJSJ5YFGhAdlR1oLCVKPoGBQiGRY9rEYdNh9tFB2FiAgAiwoNQJurC8frnRiZxqISrRQKBSYOs+CzIw2ioxARAWBRoQHYU9EKCcAoFpWoNjnLgooWF8qauJ0+EYnHokLnbeeJFiTGa3jH5Cg3JsMEtUqBzzmqQkQywKJC523nie75KdzoLbrpNSqMSTfh86MsKkQkHosKnRePz499lW2cnxIjJmVZsL20BS6vT3QUIopxLCp0Xg5UtcHrD2AkV/zEhElZFnj9AXxdwmXKRCQWiwqdlx0nWqDXKDE8mTvSxoJ0cxzSzXp8epiXf4hILBYVOi9fFTdhVJqRd9WNIRdkJ2LDoToEArzLORGJw6JC59TZ5cfOslaMz7SIjkJhdFFOIpqcXhRV2UVHIaIYxqJC57SjrAVefwDjM82io1AYFdqMMMdp8MnBetFRiCiGsajQOX1Z3ISkBC2GJcaJjkJhpFQqcEG2BR8frBMdhYhiGIsKndPmY40Yl2Hi/ikx6KLhSShr6kBxg1N0FCKKUSwqdFaN7R4crWvH+GEW0VFIgHGZZug1SnxyiKMqRCQGiwqd1VfFTQCAcRkmwUlIBK1aiQnDLPjoAIsKEYnBokJnteV4I3KS42GJ14qOQoJckpuEA9VtqGh2iY5CRDGIRYXOyB+QsOloI1f7xLjJ2YnQqZV4/0CN6ChEFINYVOiMdpe3oqXDi4tzkkRHIYH0GhUmZ1vwr30sKkQUfiwqdEbrv6lDYrwGeTaD6Cgk2LQRKThc247SRq7+IaLwYlGhfkmShPUHa3Hh8CQouSw55k3KsiBOo8L7+2tFRyGiGMOiQv06WONAjb0TF+ckio5CMqBVK3Hh8ET8k5d/iCjMWFSoX58crEOCVoUx6VyWTN2m5SWjuMGJgzVtoqMQUQxhUaF+rT9Yh0nZiVCr+CNC3SYMM8Mcp8E7u6tFRyGiGMJPITpNcUM7jtU7MYWrfegUaqUSM/KS8V5RNbr8AdFxiChGsKjQad7ZU40EXfeSVKJTXV5oRUuHF5uPNoqOQkQxgkWF+vAHJKzbU4VpI5Kh4WUf+pbhyQnISY7H23uqREchohjBTyLqY2tJE+odHlxeYBUdhWTqsgIrNh6uh93lFR2FiGIAiwr18c7uKmSY9cjnJm90BjPyUxCQgHV7OKmWiEKPRYV6OT0+rD9Yh0sLrFBwkzc6A3OcBhfnJOK17eWQJEl0HCKKckKLypYtWzBnzhxkZGRAoVDgvffeExkn5v1rXw08XQFcVpAiOgrJ3NWjU1HS2IEdZS2ioxBRlBNaVDo6OjBx4kS88MILImMQurfMX/1VGS4cnogUg050HJK5MekmZJj1eG17uegoRBTl1CLf/Prrr8f1118vMgKd9HVpM47VO/G/s0eLjkIRQKFQ4KrRqXh9RwWWzvGw3BJRyETUHBWPxwOHw9HnQcGx+ssTyEqMw9gMbplP5+fyAisUCuDNXZWioxBRFIuoorJ8+XKYzebeR1ZWluhIUaGyxYVPD9fj2rFpnERL582gV2N6Xgpe3VrOnWqJKGQiqqgsWbIEbW1tvY/KSv5LLhhWf3UCCTo1ZuRzEi0NzOzx6ahzdOLDA7WioxBRlIqooqLT6WAymfo8aGga2z1Yu70c14xJhV6jEh2HIkx2UjwmDDPj5S2lXKpMRCERUUWFgu/lLSVQKhW4fny66CgUoa4fl46DNQ4uVSaikBC66sfpdKK4uLj367KyMhQVFSEpKQnZ2dkCk8WGxnYP/r+vyzF7fDoMOqE/ChTBJg4zY1hiHF7eUoqpI5JFxyGiKCN0RGXXrl2YPHkyJk+eDABYvHgxJk+ejKVLl4qMFTN6R1PGcTSFBk+hUOCG8enYeKQBh2u5Eo+IgktoUbniiisgSdJpjzVr1oiMFROq7W68+nU5rhuXBoOeoyk0NJcWpMBm1OH5z46LjkJEUYZzVGLU0x8dRrxWhTkTMkRHoSigVirx3YkZ+OhAHY7Xt4uOQ0RRhEUlBu0ub8U/99Xi1ouyuNKHgubyQiuSDVr88bPicx9MRHSeWFRiTCAg4dfvH0ROcjxmFlhFx6EoolEpMWdiBv61vwbFDRxVIaLgYFGJMW/uqsS+yjbcdclwKJXchZaCa9ZIG1IMOjz90VHRUYgoSrCoxJCG9k488eFhzCy0YkyGWXQcikIalRK3XpSFDYfrsesE91UhoqFjUYkhj//zEJQKBe6cyj1qKHSm5yUjJzkeT354mLvVEtGQsajEiA2H6vHBgVrcdclwGPUa0XEoiikVCtw+JRt7Kuz4+GCd6DhEFOFYVGJAk9ODh9/ZjwuyLZiex51DKfTGZ5oxOcuCX//rEFxen+g4RBTBWFSinCRJWPLOAXT5A/jRZSOgUHACLYWeQqHAvOk5aHR68MLnXK5MRIPHohLl3txViQ2H63HvZSNgideKjkMxJNWkx3cnZuBPm0tR2ugUHYeIIhSLShQ7Xt+Ox/55ELNGWnFxTpLoOBSDvjsxE8kGLZasO4BAgBNriWjgWFSilNvrx4LX9sBq0GHe9BzRcShGadVK3HvpCGwva8Ffvz4hOg4RRSAWlSi19B/foKLFhQeuKoBOzW3ySZxxmWZcMyYVT390hJeAiGjAWFSi0NrtFXhrdxXmz8jFsMR40XGIcPuUbFgStPj5m0Xw+gKi4xBRBGFRiTJ7Klqx9B/f4OrRqZhZyHv5kDzoNSosvCIPB6sdWP7RYdFxiCiCsKhEkYb2Tvz0b7uRZzVg3rThouMQ9ZFvM+KHlwzH6q9O4F/7akTHIaIIwaISJTq7/PjRq7vg9QXwwFUFUKv4R0vyc82YVMzIT8ZDb+/HN9VtouMQUQTgp1kUkCQJv3x7Pw7XtON/rhmJpATul0LypFAocO+lI5Bh0WPe6h0ob+4QHYmIZI5FJQo8/1kx/rGvBvfNzEOe1SA6DtFZ6TUqPHTtKGhVStz1yg40tHeKjkREMsaiEuHe3l2FFRuO4dYLh2Ea7+NDEcIUp8GS60fB6fHh1pe+RmWLS3QkIpIpFpUI9uXxJvzynf24cpQNN0/OFB2HaECsRj0eu3EMPL4Ablm1FYdrHaIjEZEMsahEqKJKO3706i6MzzRh/owc3myQIpLN1F1WDDo1bn7hK7yxowKSxK32iejfWFQi0PH6dsz7yw5kJcXhZ1cVQq3kHyNFLku8Fo/NGYPp+Sl4eN0BLHhtD2rsbtGxiEgm+AkXYcqaOnDn/9sOc5wGv7h2FPQabo9PkU+nVuFHl43AA1cWYGtJM2b9fhOeXn8EDQ5OtCWKdQopgsdZHQ4HzGYz2traYDKZRMcJuRNNHbjt5a+hVirxyA2jYYnnMmSKPm6vH+/vr8EHB2rh80u4crQN141NwyV5yci0xImOR0RBMJDPbxaVCFHa6MTtf94G1cmSksiSQlGuw+PDVyVN2HysEaWN3futJCVokZUUh0xLHOK1aug1SkgS4A9I8PgC6Ozyw+MLwOPzw+eXoFQooFYpYIrTIDFeg3RzHEakJKAg1Yg8awLndhEJwqISZb6pbsPcv+xAvFaFX81mSaHY097ZhSO17ShvcaGxvROtri54fH50+SUoACiVgEalhFalhEalhFqlgEqhQADdJcbl9aHD40dDeyc6PH4AgDlOg4tyEnHlKBu+MzoVNpNe6O+RKJawqESRrcVN+NGru5BuicND146EUa8RHYkoYkmShPZOH040d+BYfTsO1TpwtK4dAQm4ZEQybr1wGK4fn4Z4rVp0VKKoxqISJV7bXo6l/ziIcRkmLLq6kBNniUKgvbMLeypa8cXxJhysccCoV+P2Kdm465LhyEqKFx2PKCqxqEQ4j8+PJz44jFe/Lsc1Y1Ixd1oOVEpeSycKtXpHJz49XI/PjzbA7fXjhvHpuO+KPIzNMIuORhRVWFQi2ImmDixcuwfH6ttx1yXD8Z0xaaIjEcWczi4/Nh9rxIcHatHQ7sHMQivuvzIfF+ckiY5GFBVYVCJQICDhb9vL8dRHR2CK0+CBKwuQm5IgOhZRTPMHJHxd2ox/FlWjstWNi3MSsWBWPq4otHLFENEQsKhEmIM1bVj6j4PYXd6Kq0bZcMfUbE7mI5KRgCRhT3kr/rmvBscbnBiZZsSPLxuBORMzoFVz30yigWJRiRCVLS6s2HAM7+2tRoYlDvdcmosx6ZH3+yCKFZIk4VCtAx/sr8XeSjusRh3umJKN26dkI80cO8ubAyf3rZEgQQEFtGol59GdgyR1n7OA1H3ONCoF1KrYLbkRV1ReeOEF/N///R/q6uowceJEPP/885gyZco5vy8Si4okSdhbaccrX5Tho29qYdJr8B8XDsOskTb+h04UQapaXfj4YB2+Km6G1xfAZYUp+I8LhuE7Y1IjeoVem6sLJU1OlDd3oLzZhepWN2rbOlHv6ITd3QWHuwseX+C079OrlTDo1UhO0MFq1CHDokemJR7ZyXHISU7AiBQDzPHRub2Cy+tDaWMHTpw8Z1WtbtS1uVHn6ERrRxfsLi86+zlnWpUSBp0KSQk6JBu0yLTEIcMSh+ykeAxPjscIqwEpBm1UXmaMqKLy97//HXPnzsVLL72EqVOnYuXKlXjrrbdw9OhR2Gy2s35vpBQVSZJwvMGJDYfq8c6eKpQ2diDNpMf149Mws9AKnTpy/1IjinUurw9fFTfjy+JGHKt3Qq9R4opCK64anYrLCqyyHWlp7fDiWH07jjU4cby+HcfqnTje0I5mp7f3GEu8BikGHZIStLDEaWDUq2HQqaFTq6BRK6FA92UxX0BCZ5cfHR4/2ju7YHd3oaXDi6Z2D+zurt7XSzFokW8zYmSqAQWpRhSmGlGYaoiY24E4PT4UNzhxrK69+9ydPH91bf++J5VBp4bV2H3OEuM1MOo1MOjU0GtU0J48ZwDQ5Q/A4wvA5fXD0dmFNlcXWlxeNDs9aHZ60fPBbNKrT54rAwps/z5nVqMuogtMRBWVqVOn4uKLL8Yf//hHAEAgEEBWVhb++7//Gw8//PBZv1euRcXj86OkoQP7q+zYVd6Kr0uaUW13Q6dW4uKcJFxWkIJxGWYoOYJCFFXq2jqxo6wZO8tbUdzgBADkpiTgwuGJmJRlweh0EwpTDWHbuNHp8aGq1YWKZhcqWlwoaexAaaMTxxucaOnoLiQqpQIZZj0yLHHITOy+PUG6OQ5pJj3itEP/R5Tb60edoxM1djdq7G5UnfzfGrsbgZOfPskGLfKtBuTbDMhNSUBOcgKykuIxLDEOCbrwztfr7PKj2u5GRUv3eStr6kDJyXPWU0gUAGwmPYadcs4yLHqkmeNgCEJery+Aekcnats6UW13o7rV1f2/dje6/N0nzaRXI89mQL7VgFxrAnJPnrOspHiY4+Q/chUxRcXr9SI+Ph5vv/02vve97/U+P2/ePNjtdvzjH//oc7zH44HH4+n9uq2tDdnZ2aisrAxJUZEkCV1+Cb5AAJ6uALy+AFxdPri9frR3+uDo7ILd1YVmpxeNzk5Utbq7Hy1u+CUJSgWQlRSHfJsR4zPNGJVmgjaGr0kSxRLHyW3/j9Q7cKKpAxUtrt4P5qQEDTIt8ci06JFi1CPZoIEpTgOTToN4nQp6dfe/vtVKRe8lYUkCugIB+HwSPP4A3F1+uL0+OD0+ONw+2N3dlxia2r1oaO9EfXsnnJ3+3jxatRLpZj1sRh3SzXHIMPd8uOqhVob/7yVfIIC6kx/EtW3/vrxU5/DAe8plEoNehVSjHjajHskGLZIStDDpNTDHq5GgVSPu5D2fdCoVNGpApVSi59+AktT9Pj4/4PF33wuq0+uH0+uDw92FNnfP3+EeNDo9qG/rhN3t631vtUqBVKMONqMe6SfLXLql+391qvCPhAckCQ3tnaixdxe/2rZO1DncqHd0wuX99zmL1yqRZopDqqn7nCUbtLDEa2CJ03aPiuk1SNAqEadRQ69VQadSQqdWQqM+eQuKkz93oRyxcTgcyMrKgt1uh9l89n2KhC4taWpqgt/vR2pqap/nU1NTceTIkdOOX758OR5//PHTns/KygpZxqEqB/Cl6BBEJCuVAPYJeN8SAe8ZDIcFvneZwPceiqOiA5yn9vZ2eReVgVqyZAkWL17c+3UgEEBLSwuSk5Mj+lrdqXpaZqhGiSINz0dfPB+n4znpi+ejL56PvuRyPiRJQnt7OzIyMs55rNCikpKSApVKhfr6+j7P19fXIy3t9B1ZdToddDpdn+csFksoIwpjMpn4H9UpeD764vk4Hc9JXzwfffF89CWH83GukZQeQidMaLVaXHjhhdi4cWPvc4FAABs3bsS0adMEJiMiIiI5EH7pZ/HixZg3bx4uuugiTJkyBStXrkRHRwfmz58vOhoREREJJryo3HbbbWhsbMTSpUtRV1eHSZMmYf369adNsI0VOp0Ojz322GmXuGIVz0dfPB+n4znpi+ejL56PviLxfAjfR4WIiIjoTLipBxEREckWiwoRERHJFosKERERyRaLChEREckWi4oAL7zwAnJycqDX6zF16lTs2LHjrMevXLkSI0eORFxcHLKysvDzn/8cnZ2dZ/2eSDKQ89HV1YVf//rXyMvLg16vx8SJE7F+/fowpg2tLVu2YM6cOcjIyIBCocB77713zu/ZtGkTLrjgAuh0OuTn52PNmjUhzxkuAz0ftbW1uOOOO1BYWAilUolFixaFJWe4DPR8rFu3Dt/5zndgtVphMpkwbdo0fPzxx+EJGwYDPR9ffvklZsyYgeTkZMTFxWHUqFH4wx/+EJ6wYTCYvz96fPXVV1Cr1Zg0aVLI8g0Wi0qY/f3vf8fixYvx2GOPYc+ePZg4cSKuvfZaNDQ09Hv82rVr8fDDD+Oxxx7D4cOH8corr+Dvf/87fvWrX4U5eWgM9Hw88sgj+NOf/oTnn38ehw4dwn333Yebb74Ze/fuDXPy0Ojo6MDEiRPxwgsvnNfxZWVluOGGGzBr1iwUFRVh0aJFuPfee6Pmw2ig58Pj8cBqteKRRx7BxIkTQ5wu/AZ6PrZs2YLvfOc7+PDDD7F7927MmjULc+bMidn/XhISEnD//fdjy5YtOHz4MB555BE88sgjePnll0OcNDwGej562O12zJ07F1dddVWIkg2RRGE1ZcoUaeHChb1f+/1+KSMjQ1q+fHm/xy9cuFC68sor+zy3ePFiacaMGSHNGS4DPR/p6enSH//4xz7P3XLLLdKdd94Z0pwiAJDefffdsx7z0EMPSWPHju3z3G233SZde+21IUwmxvmcj1PNnDlT+tnPfhayPKIN9Hz0GDNmjPT4448HP5Bggz0fN998s/TDH/4w+IEEG8j5uO2226RHHnlEeuyxx6SJEyeGNNdgcEQljLxeL3bv3o2rr7669zmlUomrr74aX3/9db/fM336dOzevbv3ckhpaSk+/PBDzJ49OyyZQ2kw58Pj8UCv1/d5Li4uDl9+GZv3qP7666/7nD8AuPbaa894/ii2BQIBtLe3IykpSXQUWdi7dy+2bt2KmTNnio4izOrVq1FaWorHHntMdJQzEr4zbSxpamqC3+8/bdfd1NRUHDlypN/vueOOO9DU1IRLL70UkiTB5/Phvvvui4pLP4M5H9deey1WrFiByy+/HHl5edi4cSPWrVsHv98fjsiyU1dX1+/5czgccLvdiIuLE5SM5Oj3v/89nE4nfvCDH4iOItSwYcPQ2NgIn8+HZcuW4d577xUdSYjjx4/j4YcfxhdffAG1Wr51gCMqMrdp0yY8+eSTePHFF7Fnzx6sW7cOH3zwAX7zm9+IjibEs88+i4KCAowaNQparRb3338/5s+fD6WSP8pEZ7N27Vo8/vjjePPNN2Gz2UTHEeqLL77Arl278NJLL2HlypV4/fXXRUcKO7/fjzvuuAOPP/44CgsLRcc5K/lWqCiUkpIClUqF+vr6Ps/X19cjLS2t3+959NFHcdddd/U2/vHjx6OjowM//vGP8b//+78R/QE9mPNhtVrx3nvvobOzE83NzcjIyMDDDz+MESNGhCOy7KSlpfV7/kwmE0dTqNcbb7yBe++9F2+99dZplwpjUW5uLoDuv0/r6+uxbNky3H777YJThVd7ezt27dqFvXv34v777wfQfWlQkiSo1Wp88sknuPLKKwWn7Ba5n3IRSKvV4sILL8TGjRt7nwsEAti4cSOmTZvW7/e4XK7TyohKpQIASBF+m6bBnI8eer0emZmZ8Pl8eOedd3DTTTeFOq4sTZs2rc/5A4ANGzac8/xR7Hj99dcxf/58vP7667jhhhtEx5GdQCAAj8cjOkbYmUwmHDhwAEVFRb2P++67DyNHjkRRURGmTp0qOmIvjqiE2eLFizFv3jxcdNFFmDJlClauXImOjg7Mnz8fADB37lxkZmZi+fLlAIA5c+ZgxYoVmDx5MqZOnYri4mI8+uijmDNnTm9hiWQDPR/bt29HdXU1Jk2ahOrqaixbtgyBQAAPPfSQyN9G0DidThQXF/d+XVZWhqKiIiQlJSE7OxtLlixBdXU1Xn31VQDAfffdhz/+8Y946KGHcM899+Czzz7Dm2++iQ8++EDUbyGoBno+AKCoqKj3exsbG1FUVAStVosxY8aEO37QDfR8rF27FvPmzcOzzz6LqVOnoq6uDkD3BHSz2Szk9xBMAz0fL7zwArKzszFq1CgA3cu3f//73+OBBx4Qkj/YBnI+lEolxo0b1+f7bTYb9Hr9ac8LJ3jVUUx6/vnnpezsbEmr1UpTpkyRtm3b1vtrM2fOlObNm9f7dVdXl7Rs2TIpLy9P0uv1UlZWlrRgwQKptbU1/MFDZCDnY9OmTdLo0aMlnU4nJScnS3fddZdUXV0tIHVofP755xKA0x4952DevHnSzJkzT/ueSZMmSVqtVhoxYoS0evXqsOcOlcGcj/6OHz58eNizh8JAz8fMmTPPenykG+j5eO6556SxY8dK8fHxkslkkiZPniy9+OKLkt/vF/MbCLLB/PdyKrkuT1ZIUoRfPyAiIqKoxTkqREREJFssKkRERCRbLCpEREQkWywqREREJFssKkRERCRbLCpEREQkWywqREREJFssKkRERCRbLCpEREQkWywqREREJFssKkRERCRbLCpEREQkW/8/kBc3OiUJGBIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_pa = pd.DataFrame(ATE_estimates, columns=['Estimate'])\n",
    "g = sns.kdeplot(df_pa, fill=True)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}