doc/examples/py_double_ml_ssm.ipynb

{
 "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": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================== DoubleMLData 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",
      "Selection variable: s\n",
      "No. Observations: 2000\n",
      "\n",
      "------------------ DataFrame info    ------------------\n",
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 2000 entries, 0 to 1999\n",
      "Columns: 103 entries, X1 to s\n",
      "dtypes: float64(103)\n",
      "memory usage: 1.6 MB\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from doubleml.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.DoubleMLData(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": 25,
   "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": 26,
   "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",
      "Selection variable: s\n",
      "No. Observations: 2000\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",
      "Learner ml_g_d0 RMSE: [[1.10039862]]\n",
      "Learner ml_g_d1 RMSE: [[1.11071087]]\n",
      "Learner ml_pi RMSE: [[0.42388745]]\n",
      "Learner ml_m RMSE: [[0.47222159]]\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": 27,
   "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": 28,
   "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": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABN5UlEQVR4nO3deXhU5aE/8O8smZkss2SZZLKSkJCEhBDCKqIsyiIobq31trai1rYq/lpLvbb0ehXbKnpbLa6orYLUBTekdQPZQti3EHZC9n3f15nMzPn9QUmNECBhkvfMzPfzPPM8ZHJm5stJmPlyznveVyFJkgQiIiIiGVKKDkBERETUHxYVIiIiki0WFSIiIpItFhUiIiKSLRYVIiIiki0WFSIiIpItFhUiIiKSLRYVIiIiki216ABXwul0orKyEnq9HgqFQnQcIiIiugySJKGtrQ0RERFQKi9+zMSti0plZSWio6NFxyAiIqJBKCsrQ1RU1EW3ceuiotfrAZz9ixoMBsFpiIiI6HK0trYiOjq693P8Yty6qJw73WMwGFhUiIiI3MzlDNsQOpg2NjYWCoXivNvixYtFxiIiIiKZEHpE5cCBA3A4HL1fHz9+HHPmzMEdd9whMBURERHJhdCiYjab+3z97LPPIj4+HjNmzBCUiIiIiORENmNUbDYb3n33XSxZsqTfc1ZWqxVWq7X369bW1uGKR0REAkmSBLvd3ucoPMmXSqWCWq12ydQhsikq69evR3NzM+65555+t1m+fDmeeuqp4QtFRETC2Ww2VFVVobOzU3QUGgA/Pz+Eh4dDo9Fc0fMoJEmSXJTpisybNw8ajQaff/55v9tc6IhKdHQ0WlpaeNUPEZEHcjqdyMvLg0qlgtlshkaj4QSfMidJEmw2G+rq6uBwODBq1KjzJnVrbW2F0Wi8rM9vWRxRKSkpwebNm7Fu3bqLbqfVaqHVaocpFRERiWaz2eB0OhEdHQ0/Pz/Rcegy+fr6wsfHByUlJbDZbNDpdIN+Llms9bNq1SqEhobixhtvFB2FiIhk6FLTrJP8uOpnJvwn73Q6sWrVKixatAhqtSwO8BAREZFMCG8GmzdvRmlpKe677z7RUYiIyE1UNHehqcM2bK8X6K9BpMl32F6P/kN4UZk7dy5kMp6XiIjcQEVzF65/PhPdPc5he02djxJbfjNTaFlZvXo1HnnkETQ3NwvLIILwokJERDQQTR02dPc4sXhWwrAUh4rmLry6LR9NHbYBvd4999yDd95557z7582bhw0bNlz0sbGxsXjkkUfwyCOP9N535513YsGCBZf9+oMlt0LEokJERG4p0uSLuBB/0TEu6oYbbsCqVav63DfYq1d9fX3h6+t9p5+ED6YlIiLyVFqtFhaLpc8tMDAQkiRh2bJliImJgVarRUREBH75y18CAGbOnImSkhL8+te/7l2sFzh7pMNkMvU+97JlyzBu3Di8/fbbiImJQUBAAB566CE4HA783//9HywWC0JDQ/H000/3yfTCCy8gLS0N/v7+iI6OxkMPPYT29nYAQGZmJu699160tLT0vvayZcsAnJ3L7NFHH0VkZCT8/f0xZcoUZGZmDvk+5BEVIhqQhnYr1uwpwbGKFuTXtsPo64MJIwJxXXIopieaL/0ERIRPP/0Uf/3rX7F27VqkpqaiuroaR44cAQCsW7cO6enp+PnPf46f/exnF32egoICfP3119iwYQMKCgrw/e9/H4WFhUhMTMT27duxe/du3HfffZg9ezamTJkC4Oxlwy+99BLi4uJQWFiIhx56CI899hhee+01XH311VixYgWeeOIJ5ObmAgACAgIAAA8//DBOnjyJtWvXIiIiAp999hluuOEGHDt2DKNGjRqyfcWiQkSXpcfhxKpdRXh5Sz4ckoTEMD3SIo1o7e7BxhPVWL27GLOSzFh2cypGBMv7cDzRcPniiy96P+jP+f3vfw+dTgeLxYLZs2fDx8cHMTExmDx5MgAgKCgIKpUKer0eFovlos/vdDrx9ttvQ6/XIyUlBbNmzUJubi6++uorKJVKJCUl4bnnnsO2bdt6i8q3x73ExsbiT3/6Ex544AG89tpr0Gg0MBqNUCgUfV67tLQUq1atQmlpKSIiIgAAjz76KDZs2IBVq1bhmWeeccXuuiAWFSK6pO4eBx56LxuZubWYPToM35sQBYPOp/f7kiThYEkT/rGnGPP+moXXfzIBM5NCBSYmkodZs2Zh5cqVfe4LCgpCR0cHVqxYgZEjR+KGG27AggULsHDhwgHPJxYbGwu9Xt/7dVhYGFQqVZ/J1sLCwlBbW9v79ebNm7F8+XKcPn0ara2tsNvt6O7uRmdnZ7+z/x47dgwOhwOJiYl97rdarQgODh5Q5oFiUSGii+qyOfDzfxzEvsJG/Pe8ZIyLNp23jUKhwKTYIKRFGvHK1nzc/85BvHbXeMxNvfj/Bok8nb+/PxISEs67PygoCLm5udi8eTM2bdqEhx56CH/+85+xfft2+Pj4XOCZLuy72yoUigve53SevZS7uLgYN910Ex588EE8/fTTCAoKws6dO/HTn/4UNput36LS3t4OlUqFQ4cOQaVS9fned48YuRoH0xJRvyRJwi8/OIz9RY3473lJFywp36bzUeGR2aMwMTYQD76bjW25tRfdnsib+fr6YuHChXjppZeQmZmJPXv24NixYwAAjUYDh8Ph8tc8dOgQnE4nnn/+eVx11VVITExEZWVln20u9NoZGRlwOByora1FQkJCn9ulTk9dKR5RIaJ+vbO7GJtO1eA3cxMxJtJ4WY9Rq5R4eNYovGDPxa/WHsZXv7wWUYFcTI5cr6K5S/avY7VaUV1d3ec+tVqNL774Ag6HA1OmTIGfnx/effdd+Pr6YsSIEQDOntLJysrCf/3Xf0Gr1SIkJOSK/g7nJCQkoKenBy+//DIWLlyIXbt24fXXX++zTWxsLNrb27Flyxakp6fDz88PiYmJuOuuu3D33Xfj+eefR0ZGBurq6rBlyxaMHTt2SNfqY1Ehogs6XtGCp786hRtSLZg4ImhAj1UpFXhwRgL+Z/0xPPRuNj5+cCq0atWlH0h0GQL9NdD5KPHqtvxhe02djxKB/poBP27Dhg0IDw/vc19SUhKeffZZPPvss1iyZAkcDgfS0tLw+eef9473+MMf/oBf/OIXiI+Ph9VqddkM7unp6XjhhRfw3HPPYenSpZg+fTqWL1+Ou+++u3ebq6++Gg888ADuvPNONDQ04Mknn8SyZcuwatUq/OlPf8JvfvMbVFRUICQkBFdddRVuuukml2Trj0Jy4/nrW1tbYTQa0dLSAoPBIDoOkcew2Z2Y/2IWnBLw1M2p8FEN7ixxQV07nvr8BBZNjcXjN6W4OCV5g+7ubhQVFSEuLg46na73fq71I3/9/eyAgX1+84gKEZ1nzZ5iFNV34Jnb0gZdUgAg3hyAOyZE4+1dRbg1I/KyTx8RXUqkyZfFwUtwMC0R9dHYYcOLm/NwXXKoS+ZDmZ9mQVSgL/7ns2NwON32AC4RCcKiQkR9/HXTGTglCXdMiHbJ86mVStw3bSSOlLfgg/2lLnlOIvIeLCpE1Kugrh3v7SvBbRlRMPhe/lwOl5Jk0WNmkhn/t/E0Wrt7XPa8ROT5WFSIqNfrmQUw+WkwNzXM5c99x4RodNuc+PuOIpc/N3k+N77uw2u56mfGokJEAM5eRbHucAXmj7Fc0QDa/gT5azAnJQxv7ShE4zBerUHu7dwsq52dnYKT0ECd+5kNZKbdC+FVP0QEAPhbViF8fVS4Ptn1R1POuXlcBLaersUb2wuwdMHoIXsd8hwqlQomk6l3rRo/Pz8oFArBqehiJElCZ2cnamtrYTKZzptyf6BYVIgIDe1WrN1figVjw+GrGbqJ2Qw6H8wfY8E7u4tx/7UjYdZrh+y1yHOcm6L92wvrkfyZTCaXTK/PokJEeHdvKSQA84ZhEcH5aeH4+ng11uwpxm/mJg3565H7UygUCA8PR2hoKHp6OBjbHfj4+FzxkZRzWFSIvJzd4cQH+0txdXwIDDrXXenTnwCtGjOSzFizpwQPzoyHn4ZvQ3R5VCqVyz78yH1wMC2Rl9t6uhbVrd2YPTp02F5zwRgL2rp78Omh8mF7TSJyTywqRF7uvX2lSDD7Y6Q5YNhe06zXYXJcEP62o4iz1RLRRbGoEHmx0oZOZJ2pw/Wjh+5Kn/7cNDYCpY2d2HSyZthfm4jcB4sKkRf74EAp/DQqTI0PHvbXjjcHIDEsAO/uLRn21yYi98GiQuSlHE4Jnx4qx7SEEGjVYgYoXp8chp359Siu7xDy+kQkfywqRF5qd0E9atusuHaUWViGq0YGI0Cr5mKFRNQvFhUiL/VZdgUijDrEm/2FZdColbh2VAg+OlgGq90hLAcRyReLCpEX6rTZ8fXxakxLCBE+Hfn1o8PQ1NmDDcerheYgInliUSHyQhtPVKOrx4FrR4WIjoJIky9Swg14fx9P/xDR+VhUiLzQuuwKjLboYdbrREcBAExPNGNfUSPKm7hCLhH1xaJC5GXq263YlV+PaQnij6acMzk2CFq1EusPV4iOQkQyw6JC5GU2njg7FmRSXJDgJP/hq1FhclwQPjlUDkniTLVE9B8sKkRe5qtjVUiNMA7LAoQDce0oM4obOnG4rFl0FCKSERYVIi/S2GHD3oJGTJbR0ZRzUsMNCPbXYF02Fyokov9gUSHyIhtPVEOChEmx8isqSqUC0xJC8K+cStjsTtFxiEgmWFSIvMhXx6qQEm6A0Vdep33OuTo+GK3dduzMrxMdhYhkgkWFyEs0ddiwO79Blqd9zokJ8kOkyRdfHKkSHYWIZIJFhchLbD5VA6ckz9M+5ygUClw1MggbT1Sju4dT6hMRiwqR19h8qgajwgJg8tOIjnJRU0eGoMPmQGYuT/8QkQyKSkVFBX784x8jODgYvr6+SEtLw8GDB0XHIvIo3T0ObD9ThwkxgaKjXFJkoC9GBPvhi6OVoqMQkQyoRb54U1MTpk2bhlmzZuHrr7+G2WxGXl4eAgPl/2ZK5E525deju8eJCTI+7fNtV8UF459HKtBps8NPI/RtiogEE/oO8NxzzyE6OhqrVq3qvS8uLk5gIiLPtOlkDSJMOkSafEVHuSxXjQzGhwfLkJlbhwVp4aLjEJFAQk/9/Otf/8LEiRNxxx13IDQ0FBkZGfjb3/7W7/ZWqxWtra19bkR0cU6nhE2najDeDU77nGMx6jAi2A8bj1eLjkJEggktKoWFhVi5ciVGjRqFjRs34sEHH8Qvf/lLvPPOOxfcfvny5TAajb236OjoYU5M5H4OlzWjod2GiSPc47TPORNHBGLL6VpO/kbk5YQWFafTifHjx+OZZ55BRkYGfv7zn+NnP/sZXn/99Qtuv3TpUrS0tPTeysrKhjkxkfvZfKoGRl8fjAoNEB1lQCbFBqHdasfugnrRUYhIIKFFJTw8HCkpKX3uGz16NEpLSy+4vVarhcFg6HMjoovbeqoWY6OMUCoVoqMMSEyQH8IMWmw8USM6ChEJJLSoTJs2Dbm5uX3uO3PmDEaMGCEoEZFnqWjuQm5NGzKi3Wd8yjkKhQKTYoPwzYlqOJyS6DhEJIjQovLrX/8ae/fuxTPPPIP8/Hy8//77ePPNN7F48WKRsYg8xrbTtVAqgLFRRtFRBmVSbBAaOmzILm0SHYWIBBFaVCZNmoTPPvsMH3zwAcaMGYM//vGPWLFiBe666y6RsYg8xrbTtUiy6OGvdc+5SBJCAxDo54NvTvDqHyJvJfzd66abbsJNN90kOgaRx+nucWBXQT1uy4gSHWXQlAoFxkUHYtPJGvzPjSmXfgAReRzhU+gT0dDYW9iA7h4nMqJNoqNckfEjTChu6ERhXbvoKEQkAIsKkYfadroWZr0WUYHuMRttf9IijdColNhyqlZ0FCISgEWFyENty61DepQJCoV7XZb8XVq1CmMiDdh8ipcpE3kjFhUiD1Rc34HSxk6Mc/PTPudkxATiYHETWjp7REchomHGokLkgTJza6FWKpAa4RmTIo6PCYRDkpB5hqd/iLwNiwqRB8o8U4dkix46H5XoKC4R5K/ByBB/bD7J0z9E3oZFhcjDdPc4sLegAWOjTKKjuNS4aBO2n6njLLVEXoZFhcjDHChuRLfdiXQPGZ9yTnq0Ca3dduSUcZZaIm/CokLkYbbn1iHIX4NoN78s+bsSzAHQ69TIzK0THYWIhhGLCpGHyTxTh7GRRre/LPm7lEoF0iKN2HaaA2qJvAmLCpEHqWzuQn5tu8eNTzknPcqE45WtqGuzio5CRMOERYXIg+zIq4NScXY2V090bhXorDM8/UPkLVhUiDxI1pl6xJsDEKATvt7okDD5aRBv9kdmLk//EHkLFhUiD+FwStiRV+exR1POSY8yISuvnpcpE3kJFhUiD3GsogWt3XaPHZ9yTnq0CS1dPThW0SI6ChENAxYVIg+RdaYOfhoV4kP9RUcZUvHmAPhpVBynQuQlWFSIPETWmTqkRhigVnr2P2vVv9cwYlEh8g6e/Y5G5CXauntwuLQZaZEm0VGGRVqkCYdLm9HWzdWUiTwdiwqRB9hT0ACHJPVevuvp0qOMcEgS9hQ0iI5CREOMRYXIA+zIq4fFoEOYQSc6yrAINegQbtQhK4+nf4g8HYsKkQfIOlOHMR5+WfJ3jYk0IutMvegYRDTEWFSI3FxZYydKGju95rTPOWMjjSht7ERJQ4foKEQ0hFhUiNzcjrx6KBVAaoRBdJRhlRJhgFIB7MznURUiT8aiQuTmsvLqMCpUDz+NZ06b3x8/jRqjQvXYmceiQuTJWFSI3JjDKWFXfj3GRHrX0ZRzUiMN2F3QACen0yfyWCwqRG7saHkz2rxg2vz+jIkwoqWrByerWkVHIaIhwqJC5MZ25tWfnTbfHCA6ihCjQgOgVSs5ToXIg7GoELmxrLyz0+arlArRUYRQq5QYHW7gOBUiD8aiQuSm2q32f0+b712XJX/XmAgjDhQ3orvHIToKEQ0BFhUiN7W3oAF2p+S141POGRNpgNXuRHZJk+goRDQEWFSI3NSOvDqEGbReM21+f6KD/GD09eE4FSIPxaJC5Ka2n6nDmAjvPu0DAEqFAinhBuwqYFEh8kQsKkRuqLypE8UNnV5/2uec1AgDjpW3oLW7R3QUInIxFhUiN7TTS6fN709qhBFOCThQ1Cg6ChG5GIsKkRvKyqtDfGgA/LXeNW1+f8IMWpgDtNhd0CA6ChG5GIsKkZs5O21+A9I4PqWXQqFASoQBuziglsjjsKgQuZljFS1o6erh+JTvSI0w4HR1GxraraKjEJELsagQuZkdZ+rOTpsf6i86iqykhJ8dr7O3kONUiDwJiwqRm9l+5uy0+Wol//l+W3CAFhEmHXbzMmUij8J3OiI30tbd8+9p802io8hSSjjHqRB5GhYVIjeyu6ABDknC2CgOpL2QlHAjihs6UdPaLToKEbmI0KKybNkyKBSKPrfk5GSRkYhkbUdeHcKNOq+fNr8/Kf+eV2YPL1Mm8hjCJ2FITU3F5s2be79Wq4VHIpKtrDP1GOPlqyVfjNHXB9GBvthT0IBbMyJFxyEiFxDeCtRqNSwWi+gYRLJX0tCB0sZO3DEhSnQUWRsdbuCAWiIPInyMSl5eHiIiIjBy5EjcddddKC0t7Xdbq9WK1tbWPjcib7H9TB1USgVSOdHbRaVGGFHW1IWK5i7RUYjIBYQWlSlTpmD16tXYsGEDVq5ciaKiIlx77bVoa2u74PbLly+H0WjsvUVHRw9zYiJxMnPrkGzRw1ejEh1F1pLD9QA4ToXIUwgtKvPnz8cdd9yBsWPHYt68efjqq6/Q3NyMjz766ILbL126FC0tLb23srKyYU5MJIbV7sCeggaM5fiUSzLofDAi2A97C1lUiDyB8DEq32YymZCYmIj8/PwLfl+r1UKr1Q5zKiLxDhU3oavHgbHRJtFR3ALHqRB5DuFjVL6tvb0dBQUFCA8PFx2FSFa2n6lDoJ8PRgT5iY7iFlLDDahs7kZZY6foKER0hYQWlUcffRTbt29HcXExdu/ejdtuuw0qlQo//OEPRcYikp3M3DqkRRqhUChER3ELyeEGKADs4ekfIrcntKiUl5fjhz/8IZKSkvCDH/wAwcHB2Lt3L8xms8hYRLJS09qN3Jo2pPO0z2UL0KoRG+LPcSpEHkDoGJW1a9eKfHkit7A9tw4KgBO9DdBoix57ChogSRKPRBG5MVmNUSGi823LrUVCWAAMOh/RUdxKSoQRVS3dKGvkfCpE7oxFhUjGehxO7MirR3qUSXQUt5Ns0UMB8PQPkZtjUSGSseySJrRb7RjH8SkD5q9VI47jVIjcHosKkYxty62D0dcHcSH+oqO4peRwA3YXnh2nQkTuiUWFSMa25dZibJQRSg4GHZTUcAOqW7pRyvlUiNwWiwqRTFW1dCG3ug0ZPO0zaMnheigVXPeHyJ2xqBDJ1PbcOigVQFqkSXQUt+WnOTufyr6iRtFRiGiQWFSIZGrr6VqMCtUjQCerJbnczmiLoXc+FSJyPywqRDJktTuwI68eGTEm0VHcXkq4AdWtHKdC5K5YVIhkaG9hI7p6HBgfEyg6ittLspwdp8LLlIncE4sKkQxtPVUDs16LqEBf0VHcnn/vuj8cp0LkjlhUiGRGkiRsPlWLjGgT16hxkZRwjlMhclcsKkQyk1fbjormLmTwtI/LjOY4FSK3xaJCJDNbTtVCq1YiJdwgOorHSOY4FSK3xaJCJDNbTtUgLdIIjZr/PF3FT3N23R9O/EbkfvhOSCQjjR02ZJc28bTPEBgdbsAervtD5HZYVIhkZMupGkgSMJ7zp7hcSrgBNa1WFDdwnAqRO2FRIZKRTSdrMCosACY/jegoHofzqRC5JxYVIpno7nEgK68OE3jaZ0j4adQYaQ7gOBUiN8OiQiQTu/Lr0d3jxITYINFRPNZoi57jVIjcDIsKkUxsOlmDcKMOEUad6CgeKyXCgLo2KwrrO0RHIaLLxKJCJANOp4RNp2owYUQgZ6MdQklhBo5TIXIzLCpEMnC4rAkN7TZMGMHxKUPJV6NCPMepELkVFhUiGdh4ogYmPx8khupFR/F4nE+FyL2wqBAJJkkSvj5WhQkxgVAqedpnqKVGGNDQbkN+bbvoKER0GVhUiAQ7WdWKsqYuTI7j1T7DITFMD7VSgT0cp0LkFlhUiATbeLwa/loVFyEcJjofFRJCA7A7v150FCK6DCwqRIJ9fbwa46MDoVbxn+NwSQk3YE9hI5xOjlMhkju+MxIJVFDXjrzadkziaZ9hlRphQEtXD05Xt4mOQkSXwKJCJNCG49XQqpUYG2UUHcWrJITq4aNSYHcBT/8QyR2LCpFAXx6tQkaMCVq1SnQUr6JRK5EUpueAWiI3wKJCJEhxfQdOVrXiqrhg0VG80uhwA/YVNsLucIqOQkQXwaJCJMiXx6qgVSsxLsYkOopXGhNpRLvVjuOVraKjENFFsKgQCcLTPmKNNPvD10fFcSpEMseiQiQAT/uIp1YqkWzRYxfnUyGSNRYVIgF42kceUiOMOFjchO4eh+goRNQPFhUiAT4/UsnTPjKQGmmA1e5EdmmT6ChE1A8WFaJhll/bhtPVbbg6PkR0FK8XE+QHg06N3fm8TJlIrlhUiIbZv3Iq4a9RYVy0SXQUr6dUKJASYcAuDqglki0WFaJhJEkS/plTiYmxQfDh2j6ykBphxNGyFrR194iOQkQXwHdKomF0rKIFJY2duDqeV/vIRWqEAQ5Jwv6iRtFRiOgCZFNUnn32WSgUCjzyyCOioxANmX/lVMLo64PUCK7tIxcWgw7mAC128jJlIlmSRVE5cOAA3njjDYwdO1Z0FKIh43RK+PxoJSbHBUGlVIiOQ/+mUCgwJtKAHXksKkRyJLyotLe346677sLf/vY3BAYGio5DNGT2FjWgptWKabzaR3bGRBqRX9uOmtZu0VGI6DsGVVQKCwtdFmDx4sW48cYbMXv27Etua7Va0dra2udG5C7WH65AmEGLxLAA0VHoO8b8+1TcTh5VIZKdQRWVhIQEzJo1C++++y66uwf/P5C1a9ciOzsby5cvv6ztly9fDqPR2HuLjo4e9GsTDafuHge+OlaNaQkhUCh42kduDL4+iAvx5zgVIhkaVFHJzs7G2LFjsWTJElgsFvziF7/A/v37B/QcZWVl+NWvfoX33nsPOp3ush6zdOlStLS09N7KysoGE59o2G0+VYN2qx3X8LSPbKVGGLAjrw6SJImOQkTfMqiiMm7cOLz44ouorKzE22+/jaqqKlxzzTUYM2YMXnjhBdTV1V3yOQ4dOoTa2lqMHz8earUaarUa27dvx0svvQS1Wg2H4/y1N7RaLQwGQ58bkTv4LLsCCaEBCDf5io5C/UiLNKK+3YYzNe2ioxDRt1zRYFq1Wo3bb78dH3/8MZ577jnk5+fj0UcfRXR0NO6++25UVVX1+9jrr78ex44dQ05OTu9t4sSJuOuuu5CTkwOVimugkGdo7LBh+5k6DqKVuWSLAT4qBXbkXfo/WkQ0fK6oqBw8eBAPPfQQwsPD8cILL+DRRx9FQUEBNm3ahMrKStxyyy39Plav12PMmDF9bv7+/ggODsaYMWOuJBaRrPwrpwISwEneZE6jViLZwsuUieRGPZgHvfDCC1i1ahVyc3OxYMECrFmzBgsWLIBSebb3xMXFYfXq1YiNjXVlViK39MmhcmREm2Dw9REdhS5hbJQRnx4qR3ePAzofHtUlkoNBFZWVK1fivvvuwz333IPw8PALbhMaGoq33nprQM+bmZk5mDhEspVb3Ybjla1YMidRdBS6DGOjTHhvXykOFDfi2lFm0XGICIMsKps2bUJMTEzvEZRzJElCWVkZYmJioNFosGjRIpeEJHJXn2aXQ69TI4MrJbuF6EBfBPlrkHWmjkWFSCYGNUYlPj4e9fXnn8dtbGxEXFzcFYci8gR2hxOfZpdjWnwI1Fwp2S0oFAqkRRqx/QwH1BLJxaDePfubZ6C9vf2y50Qh8nQ78urR0G7D9ET+z9ydpEcZcaamHVUtXaKjEBEGeOpnyZIlAM7+r+OJJ56An59f7/ccDgf27duHcePGuTQgkbv68GAZYoL8EBvsd+mNSTbGRBqhALDjTD1+MImzXxOJNqCicvjwYQBnj6gcO3YMGo2m93sajQbp6el49NFHXZuQyA01tFux+WQNfjQlhlPmuxm9zgfxoQHYfqaORYVIBgZUVLZt2wYAuPfee/Hiiy9yZliifnx2uAIKBXBNAid5c0djo4zYfLIGdoeT44uIBBvUv8BVq1axpBD1Q5IkfHigDBNGBEKv49wp7mhclAmt3XbklDWLjkLk9S77iMrtt9+O1atXw2Aw4Pbbb7/otuvWrbviYETu6kh5C/Jq2/G98cmio9AgxZsDYNCpsS23FhNjg0THIfJql11UjEZj77l2o9E4ZIGI3N2HB8oQEqBBWiT/nbgrpVKBsVEmbD1di/+ex8JJJNJlF5VVq1Zd8M9E9B8dVjv+lVOBeakWKJUcROvOxkWb8Mq2fFS3dMNi5LQLRKIMaoxKV1cXOjs7e78uKSnBihUr8M0337gsGJE7+uJoJTptDsxMChUdha5QepQJSgWQmVsrOgqRVxtUUbnllluwZs0aAEBzczMmT56M559/HrfccgtWrlzp0oBE7uT9/aVIjzbCrNeKjkJXKECnxqgwPbaxqBAJNaiikp2djWuvvRYA8Mknn8BisaCkpARr1qzBSy+95NKARO7iVFUrjpS14LqkMNFRyEXGRZmwI68eNrtTdBQirzWootLZ2Qm9Xg8A+Oabb3D77bdDqVTiqquuQklJiUsDErmLtftLYfL1QcYIk+go5CIZMSZ02hzYX9QoOgqR1xpUUUlISMD69etRVlaGjRs3Yu7cuQCA2tpazq9CXqnL5sC6wxWYnmiGWskJwjxFTJAfQgI02HyqRnQUIq81qHfUJ554Ao8++ihiY2MxZcoUTJ06FcDZoysZGRkuDUjkDr44Wom2bjuuS+YgWk+iUCgwPiYQm07W9LsYKxENrQFNoX/O97//fVxzzTWoqqpCenp67/3XX389brvtNpeFI3IX7+4rQXqUEWEGXsbqaSaMCMQ3J2uQW9OGZAuPGBMNt0EVFQCwWCywWCx97ps8efIVByJyNycqW3CkrAVL5iSKjkJDYHS4Ab4+Kmw+WcOiQiTAoIpKR0cHnn32WWzZsgW1tbVwOvuOiC8sLHRJOCJ38P6+UgT5azA+JlB0FBoCPiolxkYZ8c3JGjx83SjRcYi8zqCKyv3334/t27fjJz/5CcLDw7mMPXmtdqsdnx2uwA1jLFBxJlqPNWFEIF7LLEBtazdCeXqPaFgNqqh8/fXX+PLLLzFt2jRX5yFyK+sPV6C7x4HrOBOtRxsXfXaW2s2navGjKTGi4xB5lUFd9RMYGIigIK4oSt5NkiT8Y08JJowIRHAAZ6L1ZHqdD5ItBmw8US06CpHXGVRR+eMf/4gnnniiz3o/RN7mYEkTcmvaMHs0Z6L1BpNiA7Ervx6t3T2ioxB5lUGd+nn++edRUFCAsLAwxMbGwsfHp8/3s7OzXRKOSM7+sacY4UYdxkQaRUehYTApNgjv7CnBttO1uGVcpOg4RF5jUEXl1ltvdXEMIvdS12bFV8eq8cPJMVByMLlXCA7QIt7sjw3Hq1lUiIbRoIrKk08+6eocRG7lo4NlUCkVmJ5oFh2FhtHE2CD8M+fsAGqdj0p0HCKvMOhFSZqbm/H3v/8dS5cuRWPj2QW7srOzUVFR4bJwRHJkdzjx7t4STB0ZjADtoOdMJDc0OTYI3T1OZJ2pEx2FyGsMqqgcPXoUiYmJeO655/CXv/wFzc3NAIB169Zh6dKlrsxHJDubT9WiqqUbc1Mtl96YPEqEyRdRgb7YwKt/iIbNoIrKkiVLcM899yAvLw863X8mP1qwYAGysrJcFo5Ijt7ZXYykMD3iQvxFRyEBJscGYdOJGtjszktvTERXbFBF5cCBA/jFL35x3v2RkZGorub/NMhz5dW0YU9hA+ak8JJkb3XVyGC0We3Ymc/TP0TDYVBFRavVorW19bz7z5w5A7OZgwvJc63ZUwKTrw+mxHHCQ28VFXj29M8XR6tERyHyCoMqKjfffDP+8Ic/oKfn7MRHCoUCpaWl+O1vf4vvfe97Lg1IJBet3T345FA5rhsdCrVq0OPQyc0pFApMiQvCNydqYLU7RMch8niDerd9/vnn0d7eDrPZjK6uLsyYMQMJCQnQ6/V4+umnXZ2RSBY+PlgOm8PJmWgJU+KC0W61Y8eZetFRiDzeoK6tNBqN2LRpE3bt2oUjR46gvb0d48ePx+zZs12dj0gWnE4J7+wuxpS4IAT6aUTHIcGig/wQHeiLz49WYjbHKxENqQEXFafTidWrV2PdunUoLi6GQqFAXFwcLBYLJEmCgrN0kgfKPFOL0sZO3H9NnOgoJBNTRgbjq2NVnPyNaIgN6NSPJEm4+eabcf/996OiogJpaWlITU1FSUkJ7rnnHtx2221DlZNIqFW7ihFv9kdCaIDoKCQTV8cHo9PmwNbTtaKjEHm0AR1RWb16NbKysrBlyxbMmjWrz/e2bt2KW2+9FWvWrMHdd9/t0pBEIuXXtmNHXj0enBHPI4bUK9zoiwSzP9YfrsCCtHDRcYg81oCOqHzwwQf4/e9/f15JAYDrrrsOv/vd7/Dee++5LByRHKzeXQSTrw+mxgeLjkIyMzU+BNtya9HS2SM6CpHHGlBROXr0KG644YZ+vz9//nwcOXLkikMRyUVL59lLkq8fHQYfXpJM3zE1PhgOp4Svj3NOFaKhMqB33sbGRoSF9T/CPSwsDE1NTVccikgu1h4ohcMpYfboUNFRSIYC/TRIjTBifU6l6ChEHmtARcXhcECt7n9Yi0qlgt1uv+JQRHJgdzjxzu5iTB0ZDBMvSaZ+TEsIxr7CBlS1dImOQuSRBjSYVpIk3HPPPdBqtRf8vtVqHdCLr1y5EitXrkRxcTEAIDU1FU888QTmz58/oOchGgqbTtagsqUbD183SnQUkrFJsUF4e2cx1h+uxIMz40XHIfI4AyoqixYtuuQ2A7niJyoqCs8++yxGjRoFSZLwzjvv4JZbbsHhw4eRmpo6kGhELvf3nUUYHc5Vkuni/DRqTIwNxCeHyvDAjJG8MozIxRSSJEmiQ3xbUFAQ/vznP+OnP/3pJbdtbW2F0WhES0sLDAbDMKQjb3GkrBm3vLoLS+YkYlIsFyCkiztS1oxnN5zGPxdPQ3q0SXQcItkbyOf3oKbQHwoOhwMff/wxOjo6MHXq1AtuY7Va+5xeutAKzkSu8NbOQoQZtJgQEyg6CrmBtEgjgvw1+DS7nEWFyMWEX2957NgxBAQEQKvV4oEHHsBnn32GlJSUC267fPlyGI3G3lt0dPQwpyVvUNnchS+PVeOGVAuUSh7Gp0tTKhWYFh+Mf+ZUckVlIhcTXlSSkpKQk5ODffv24cEHH8SiRYtw8uTJC267dOlStLS09N7KysqGOS15g3f2FEOrVmJGIi9Jpss3PdGMlq4ebOOU+kQuJfzUj0ajQUJCAgBgwoQJOHDgAF588UW88cYb522r1Wr7veKIyBU6rHa8v68Us5JC4avhQnN0+aIC/ZBg9seHB8pwwxhOqU/kKsKPqHyX0+kc8GXORK7y0cEydFjtuGGMRXQUckPTE0Ox/Uwdqlu6RUch8hhCi8rSpUuRlZWF4uJiHDt2DEuXLkVmZibuuusukbHISzmcEt7aWYSrRgYjJIBH7mjgpiUEw0elxKfZ5aKjEHkMoUWltrYWd999N5KSknD99dfjwIED2LhxI+bMmSMyFnmpb05Uo7ypiyvh0qD5adSYHBeEDw+UwemU1cwPRG5L6BiVt956S+TLE/Xx5o5CjA7XI94cIDoKubFZSaH4wxcnsa+okStuE7mA7MaoEIlwqKQRh0ubeTSFrliyRY9wow5rD5SKjkLkEVhUiAC8mVWISJMvxnOCN7pCCoUCMxPN+PpYNZo7baLjELk9FhXyekX1HfjmRA0WpIVDyXVayAVmJIXCKUn4NLtCdBQit8eiQl7vrZ2FMPj64JqEENFRyEMYfX0wMTYQ7+0rgcyWUyNyOywq5NUa2q34+GA55qaEQaPmPwdyneuTw1BY14H9RY2ioxC5Nb4zk1d7Z08JAGBOSpjgJORpUiMMCDfq8P4+DqoluhIsKuS1Om12vLO7GLOSQqHX+YiOQx5GoVDguuRQfHW8Co0dHFRLNFgsKuS1PjpQhrbuHixI43T5NDSmJ5oBAB8f5AKqRIPFokJeye5w4s0dhbg6PgRmvU50HPJQBp0Ppo4Mxj/2lnCmWqJBYlEhr/TlsSpUNnfjprGc4I2G1pyUMJQ3dWF7Xp3oKERuiUWFvI4kSXgtswDpUUaMCPYXHYc8XLw5AHEh/vjHvwduE9HAsKiQ18k8U4fc6jbcnB4hOgp5AYVCgdmjw7DtdC3KGjtFxyFyOywq5HVWbivAqNAAjA43iI5CXuLq+GD4aVR4j5cqEw0Yiwp5lUMljdhf3IiF6RFQcLp8GiY6HxVmJJqx9kApunscouMQuRUWFfIqKzMLEGnyxYQRXHyQhtecFAuaO3vw+ZFK0VGI3AqLCnmN3Oo2bD5Vi4XpXHyQhp/FqMO4aBPe2V3M9X+IBoBFhbzGysx8hARoMC2eiw+SGHNTwnC8shU5Zc2ioxC5DRYV8gpljZ34/EgVbkwLh1rFX3sSIz3ahDCDFqt3F4uOQuQ2+I5NXuGNrAL4a1WYlRwqOgp5MaVCgTmjLfjyaBVq27pFxyFyCywq5PFq27rx0YFyzEu1QKtWiY5DXm5mkhkqpQIf7OP6P0SXg0WFPN5bO4qgVikwL5WLD5J4/lo1rkkIwbv7SmCzO0XHIZI9FhXyaM2dNvxjbwnmpITBX6sWHYcIADAv1YK6Nis2nKgWHYVI9lhUyKOt3l0Mh1PC/DFcfJDkIzrID6kRBqzeVSQ6CpHssaiQx2q32vH2riLMSgqF0ddHdByiPualWJBd2ozjFS2ioxDJGosKeaz39pag0+rATWN5NIXkZ/yIQJgDNFi9m0dViC6GRYU8UnePA2/uKMT0RDOCA7Si4xCdR6U8u6ryv3Kq0NBuFR2HSLZYVMgjfXigDE0dNtycHiE6ClG/zs3rs/YAL1Um6g+LCnkcm92JlZkFuDo+BGEGneg4RP3S63xwdXww/rGnBHYHL1UmuhAWFfI4nx0uR3VrN24Zx6MpJH/zxlhQ3dqNTSdrREchkiUWFfIodocTr24rwOTYIEQF+omOQ3RJscH+SLbouf4PUT9YVMijfHG0CqWNnbg1I1J0FKLLNjclDPuKGnG6ulV0FCLZYVEhj+F0Snh5ax7Gx5gQF+IvOg7RZZsUF4RAPx+s2VMiOgqR7LCokMfYcKIaBXUduHUcj6aQe1Erlbh+dBjWZZejpbNHdBwiWWFRIY8gSWePpqRFGjEqTC86DtGAXZ8cCrtDwseHeKky0bexqJBH2Hq6Fqeq2jg2hdyWyU+DyXFB+MeeEjidkug4RLLBokJuT5IkvLQ1D8kWPUZbeDSF3NfcFAtKGjuxI79edBQi2WBRIbe3K78BR8pacOu4SCgUCtFxiAYtMSwAscF+WMNLlYl6saiQ23t5ax7izf4YG2UUHYXoiigUCsxOCcPW07Uoa+wUHYdIFlhUyK3tL2rEvqJG3MKjKeQhpsWHwE+jwnv7SkVHIZIFFhVya69szUNMkB8mjAgUHYXIJXQ+KlybaMaHB0phtTtExyESjkWF3NaRsmZk5dXj1nERUPJoCnmQ2aPD0NTZgw3Hq0VHIRJOaFFZvnw5Jk2aBL1ej9DQUNx6663Izc0VGYncyCtb8xFh1GFKXLDoKEQuFWnyRWqEgTPVEkFwUdm+fTsWL16MvXv3YtOmTejp6cHcuXPR0dEhMha5gdPVrdh0qgY3j4uEUsmjKeR55owOw6GSJq7/Q15PLfLFN2zY0Ofr1atXIzQ0FIcOHcL06dMFpSJ38OrWfJj1WkxL4NEU8kwTYgMR6OeDd/eW4E+3pomOQySMrMaotLS0AACCgoIu+H2r1YrW1tY+N/I+hXXt+OJoFRaODYdaKatfYSKXUSuVmJUUinXZFWi32kXHIRJGNu/yTqcTjzzyCKZNm4YxY8ZccJvly5fDaDT23qKjo4c5JcnByswCmPx8MCMxVHQUoiE1KzkU3T0O/CunUnQUImFkU1QWL16M48ePY+3atf1us3TpUrS0tPTeysq4eJe3KW/qxLrDFbgxLQIatWx+fYmGREiAFuOiTXhvHwfVkveSxTv9ww8/jC+++ALbtm1DVFRUv9tptVoYDIY+N/Iub2YVwk+jwvWjeTSFvMP1yWE4UdmKo+XNoqMQCSG0qEiShIcffhifffYZtm7diri4OJFxSOZq27qxdn8Zbki1QOejEh2HaFiMizYhJECD9zlTLXkpoUVl8eLFePfdd/H+++9Dr9ejuroa1dXV6OrqEhmLZOrtncVQKRWYm2oRHYVo2CiVCsxKCsX6nAq0dveIjkM07IQWlZUrV6KlpQUzZ85EeHh47+3DDz8UGYtkqKWzB//YU4w5KWEI0Aq9qp5o2M1MCoXN7uSgWvJKQt/xJUkS+fLkRtbsKUaPQ8L8MTyaQt4nyF+D8TGBeH9fKX581QjRcYiGlSwG0xJdTKfNjrd2FWFGkhkmP43oOERCzEoOxcmqVhwrbxEdhWhYsaiQ7K3dX4bWrh4sHBsuOgqRMOOiTAgO0OD9/RxUS96FRYVkzWZ34s2sQkxLCIFZrxMdh0gYpVKBGYlm/DOnAh2cqZa8CIsKydr6nApUt3Zj4dgI0VGIhJuVFIoumwOfH+GgWvIeLCokW06nhJWZBZg4IhDRQX6i4xAJFxKgRXq0EWsPcFZu8h4sKiRb35ysRlF9B25O59EUonNmJoUip6wZudVtoqMQDQsWFZIlSZLw2rYCpEYYMCpMLzoOkWxMiAmE0dcHH/KoCnkJFhWSpT0FDTha0cKxKUTfoVYpcU1CCNYdLofV7hAdh2jIsaiQLL2WWYDYYD+MjTKKjkIkO7OSQtHc2YNNJ2tERyEaciwqJDvHK1qwM78eC9MjoFAoRMchkp3IQF8kWfRYu5+nf8jzsaiQ7KzMLECYQYspccGioxDJ1sxEM3bl16O8qVN0FKIhxaJCslLS0IGvj1dhQVo4VEoeTSHqz1Ujg6HzUeGTQ+WioxANKRYVkpU3swqh1/lgZmKo6ChEsqbzUeGqkUH46GAZnE4u8Eqei0WFZKOuzYqPD5ZjbkoYNGr+ahJdysykUFQ2d2N3QYPoKERDhp8GJBurdxdBqQTmplhERyFyC6NCAxBp8sWHB7hQIXkuFhWShXarHWv2lOC6pFAE6NSi4xC5BYVCgZlJZmw8UYPmTpvoOERDgkWFZGHt/lJ02hxYkBYuOgqRW7kmIQR2pxP/4kKF5KFYVEg4m92Jv+8owrT4YAQHaEXHIXIrJj8NxscEckp98lgsKiTc50cqUd3ajZs4XT7RoMxMCsWJylacqGwRHYXI5VhUSCinU8Lr2wswPsaE6CA/0XGI3NK4aBMC/Xzw8UHOqUKeh0WFhMo8U4u82nYeTSG6Aiqlonehwu4eLlRInoVFhYR6PbMQo8ICkGzRi45C5NZmJIWitcvOhQrJ47CokDDZpU3YX9yIhWlcfJDoSkWafJEUpsfHBzmoljwLiwoJ88b2AkQYdZgwIlB0FCKPMCPRjB159aho7hIdhchlWFRIiIK6dnxzogYLxoZDycUHiVziqpHB0Poo8SkXKiQPwqJCQvwtqxAmPx9cm2AWHYXIY/hqVJgSF8yFCsmjsKjQsKtt7can2eW4IdXCxQeJXGxmkhnlTV3YW8iFCskz8FOCht1bu4rgo1JidkqY6ChEHicpTI8Io44z1ZLHYFGhYdXa3YP39pbi+uRQ+Gm4+CCRqykUCsxICsXXx6vR0tkjOg7RFWNRoWH13t5SWO0OzOfig0RDZvqoEDgkCetzKkRHIbpiLCo0bLp7HPj7zkJcO8qMQD+N6DhEHuvsQoUmfLC/FJLEQbXk3lhUaNh8cqgcTR02LOR0+URDbmZSKE5Xt+F4RavoKERXhEWFhoXd4cQb2wswOS4IFqNOdBwij5ceZUKQvwZrD5SKjkJ0RVhUaFh8dbwaZU1duDk9UnQUIq+gUiowfZQZ/8ypRKfNLjoO0aCxqNCQkyQJr23Lx9hII+JC/EXHIfIas5LM6LDa8cXRKtFRiAaNRYWG3NbTtThd3YZbMng0hWg4hRp0SIsy4oP9PP1D7otFhYaUJEl4eWs+ksL0GG3Ri45D5HWuSw7F4dJm5Fa3iY5CNCgsKjSk9hY2IqesGbeMi4BCwcUHiYbbhBGBMPn58KgKuS0WFRpSr27LR2ywH8ZFm0RHIfJKaqUS00eZ8Wl2Obp7HKLjEA0YiwoNmcOlTdiZX4+b0yN5NIVIoOuSQ9HWzUG15J6EFpWsrCwsXLgQERFnTwusX79eZBxysZe25CEq0BdT4oJERyHyamEGHdKjjPjH3mLRUYgGTGhR6ejoQHp6Ol599VWRMWgIHK9owbbcOtwyLhJKJY+mEIk2e3QYjpS14HhFi+goRAMidPna+fPnY/78+SIj0BB5aUseLEYdpo4MFh2FiABkxAQiOECD9/aVYPntY0XHIbpsbjVGxWq1orW1tc+N5OdUVSu+OVmDW8dFQMWjKUSyoFIqcF1SKD47XIHW7h7RcYgum1sVleXLl8NoNPbeoqOjRUeiC3hx8xmEGbSYlhAiOgoRfcus5FD0OCR8eqhcdBSiy+ZWRWXp0qVoaWnpvZWVlYmORN9xorIFG07U4LaMSKiVbvXrReTxAv00mBIXhNW7i+F0SqLjEF0Wt/ok0Wq1MBgMfW4kLys2nR2bck2CWXQUIrqAeakWlDR0YntenegoRJfFrYoKydux8hZsOlWD2zMiOTaFSKZGhQYg3uyP1buKRUchuixCi0p7eztycnKQk5MDACgqKkJOTg5KSznVszv6yze5iDD54up4jk0hkiuFQoG5KRZsP1OHgrp20XGILkloUTl48CAyMjKQkZEBAFiyZAkyMjLwxBNPiIxFg7C/qBHbz9ThjglRPJpCJHNT44Nh9PXhURVyC0LnUZk5cyYkiQO63J0kSXhuw2mMDPHHZM5CSyR7PiolZo8Ow8cHy7BkTiIC/TWiIxH1i2NU6Ipl5tbhUEkT7pgYDSXX9CFyC3NTw+CUgHf3loiOQnRRLCp0RRzOs0dTRofrkR5lFB2HiC6TQeeD6YlmrN5dzFWVSdZYVOiKrMsux+nqNvxocgxXSCZyMzemhaOxw4bPDleIjkLULxYVGrQumwN/+SYXV40MQkKoXnQcIhogi1GHSbFBeCOrAA5OAEcyxaJCg/b2riI0tNvwX5NiREchokG6eVwEius78eWxKtFRiC6IRYUGpbatG69ty8fslDCEGXSi4xDRIMWbAzAu2oiXt+RxWn2SJRYVGpQ/b8iFUqHA9zKiREchoit0W0YU8mrb8c3JatFRiM7DokIDdqSsGR8fKsf3J0YhQCd0Kh4icoHEMD1SIwx4cUse57Yi2WFRoQFxOiU8+a8TiAnyw/XJYaLjEJGL3J4RiVNVbdhwnEdVSF5YVGhAPs0uR05ZMxZNHcGp8ok8SEqEEWOjjPi/jbmwO5yi4xD1YlGhy9bUYcPTX53CNQkhSIng5G5Enua/JsWgqL4DnxwqFx2FqBeLCl225V+fQo/dibum8HJkIk8UF+KPq+OD8dfNZzhbLckGiwpdlgPFjfjoYDnunBQDkx8XMCPyVHdMiEZ9uw1/31EoOgoRABYVugzdPQ489slRjAoLwPWjQ0XHIaIhZDHqcEOqBa9sy0dlc5foOEQsKnRpf910BuVNnfjFtfFcHZnIC9w+PhK+Pio8/eUp0VGIWFTo4g6XNuFvOwrx/fFRiAz0FR2HiIaBn0aNH06OwZfHqrA7v150HPJyLCrUry6bA7/5+AjiQvxx49gI0XGIaBhdkxCCJIsev//sGLpsHFhL4rCoUL+e/uokKpq68OCMBM6ZQuRlFAoFfn7tSFQ2d+PPG3NFxyEvxqJCF7TpZA3e3VuKu6aM4CkfIi8VYfLFDyZGY9WuIuwvahQdh7wUiwqdp6qlC499cgQTRwRiNq/yIfJq88dYkGjR4zcf5aClq0d0HPJCLCrUh83uxEPvZUOpUOBn00dCwat8iLyaUqnAgzPi0dhpw28+yoHTyUULaXixqFAfy78+hWPlLfjV9aNg0PmIjkNEMhBm0OGhGQnYfKoWr2cViI5DXoZFhXqtP1yBVbuK8eOrRmBUmF50HCKSkfEjAnHruEj8ZWMuNp2sER2HvAiLCgEADpU04bFPjmJ6YgjmpoSJjkNEMnTHhChMjA3C4veysa+wQXQc8hIsKoTypk78fM1BxJv9cf81HJdCRBemVCrw8KwEJIYF4KfvHMSx8hbRkcgLsKh4ucYOG+5+ez/UKgUemZMIHxV/JYiofz4qJZbMSUK4UYc739yDHXl1oiORh+Onkhdrt9pxz9v70dBuw29vSObgWSK6LL4aFX6/YDQSw/S4Z9UBfHSwTHQk8mAsKl6qy+bAz945iPy6dvz2hmSEGzmpGxFdPp2PCr+Zm4gZiWY89slR/OqDw5xnhYYEi4oX6rTZcd/qA8gubcKjc5MQF+IvOhIRuSG1UomfXTsSD89KwOZTNZj31yx8cqgcDs61Qi6kkCTJbX+jWltbYTQa0dLSAoPBIDqOW2i32vHT1QdwpLwZv52XjORw7jciunJ1bVa8t68E+4oakRgWgHunxWFhegQCtGqXvUaXzYG6NivqO6xo6rCh3WpHh9WBHocTDqcEpQLQqFXw1Shh8tXA5OcDi1GHUL2O65XJzEA+v1lUvEhtWzfuXXUARfUd+O95SUi2cJ8RkWvl17ZhXXYFjpQ3Q6NWYvooM6bGByMjJhCxwX4w+Wku+Dib3YmGDiuqWrpR3dKN8qZOVDR1oby5C2WNnahu6UZrt/28xykAqFUKKBUKSBLQ43Diux9qaqUCESZfJIQGIDFMjzGRBqRFGhET5MerHAVhUaHzFNS1Y9Hb+9FhteO3NyRjRDBP9xDR0Glot2JHXj2OVjQjr7YddsfZjxo/jQoBWjV0PipIkgSbw4m2bjs6bY4+j/f1USFEr0GIvxbBAVoEB2gQ7K+ByU8Dg04Nvc4HfhoVtGpln7Jx7jk7rA60dfegocOGhnYrqlutqGzuQllTJxrabQCAIH8NJowIxFUjg3F1fDCSwvRQ8sjLsGBRoT62nKrBr9bmwOTng8fmJcOs14qORERexGZ3oqK5CzWt3ahvt6K7xwGr3QmlQgG1SgFfHxX0Oh8YdGoE+WsQ7K+Fv1Y1ZEc7Wrt6UFDXjrzadpypacOZmjb0OCQEB2gwI9GMmUmhmJFohtGXV0IOFRYVAgA4nBJe3pqHFzfnYcKIQDw4Mx5+GtedLyYi8gQ2uxNnatpwtLwZR8pbUNrYCbVSgclxQZiXasHc1DBeGeliLCqEiuYuPLL2MA6VNOH28VG4LSMSSp6LJSK6pPp2K7JLm3CopAknKlvhcEoYF23CjWnhmJ9mQVSgn+iIbo9FxYtJkoQPD5Th6S9PQeujxOKZCbyyh4hokDqsdmSXNmF/USOOlDejxyEhPdqIhWMjMD8tHJEmHmkZDBYVL5VX04Yn/nkCewobMCPRjB9fNcKllwYSEXmzLpsD2aVN2FvY0FtaMmLOHmlZkBaOCJaWy8ai4mUaO2x4aUse/rGnBCF6De6bFoexUSbRsYiIPFanzY7s0mbsLWjA0Yq+pYVHWi6NRcVLNHXY8PedhVi1qxiSBNyaEYn5YyxcWJCIaBh12uw4VNL39FBapBE3jLFgbkoYEkIDOF/Ld7CoeLjCunas2lWMjw+VQZKAeakW3Dg2nIsKEhEJ1mmzI6esGQeKG5FT1ozuHidGBPlhdkoYrksOxaTYIGjU/M8ki4oH6rDa8c3JaqzdX4Z9RY0w+vpg9ugwzE0Jg4HX+hMRyY7N7sTxyhYcKmlCTlkzGjts8PVR4eqEYEwfZcbV8cFee7TF7YrKq6++ij//+c+orq5Geno6Xn75ZUyePPmSj/P0otLQbkVWXh02nqjBttO1sNqdSI0wYEaiGVPigtnKiYjchCRJKGnsxJGyZhyraEFudRvsTgkhARpMiQvG5LggTBgRiGSLHmovOH3vVkXlww8/xN13343XX38dU6ZMwYoVK/Dxxx8jNzcXoaGhF32spxWVmtbus4cMixqxt7ABJypbIQFIMPtjclwwrhoZBLNeJzomERFdoe4eB87UtOFEZStyq9tQUNcOu1OCzkeJ1Agj0iKNGBNpRLJFj4TQAOh8VKIju5RbFZUpU6Zg0qRJeOWVVwAATqcT0dHR+H//7//hd7/73UUf645FxemUUNduRUlDJ4rrO1BQ147c6jacrGpFbZsVAGAO0CLZokdqpAHpUaZ+F/EiIiLPYLM7UVTfgbzaNhTWd6C4vgNVLd0AAKUCiAr0Q0JoAGKD/TEi2A/RQb6INPkhwqRDgFbtdqePBvL5LXSSDZvNhkOHDmHp0qW99ymVSsyePRt79uw5b3ur1Qqr1dr7dUtLC4Czf+GhIkkSnBJgd55dRrzHIcHhcML271uP3YmuHge6bU502h3otNrR3m1Hm7UHrV12NHX2oLnThto2K+r+fbM7/9MNzQEaRJh8MT7cF3FjghEX4ocg/2+txeOwoqPNeoFkRETkSaL8gag4PWbF6QEAXT0OVDR3ory5C9XN3aiqb8Lx4irUt9n6fI5ofZQwB2gR7K9BUIAGgb4aGP18oNepodep4a9Rw1ejhr/27CKOOrUKOh8lfNRKaFRK+KiUUKuUUCvPrr2kUiigUiqGtPyc+9y+nGMlQotKfX09HA4HwsLC+twfFhaG06dPn7f98uXL8dRTT513f3R09JBlHGplALJFhyAiIreWLzrAILW1tcFoNF50G7eatnTp0qVYsmRJ79dOpxONjY0IDg5GW1sboqOjUVZW5jangTxJa2sr978g3PficN+Lw30v1pXuf0mS0NbWhoiIiEtuK7SohISEQKVSoaamps/9NTU1sFgs522v1Wqh1Wr73GcymQCg9xCVwWDgL61A3P/icN+Lw30vDve9WFey/y91JOUcoddAaTQaTJgwAVu2bOm9z+l0YsuWLZg6darAZERERCQHwk/9LFmyBIsWLcLEiRMxefJkrFixAh0dHbj33ntFRyMiIiLBhBeVO++8E3V1dXjiiSdQXV2NcePGYcOGDecNsL0UrVaLJ5988rxTQzQ8uP/F4b4Xh/teHO57sYZz/wufR4WIiIioP54/Ty8RERG5LRYVIiIiki0WFSIiIpItFhUiIiKSLbcqKq+++ipiY2Oh0+kwZcoU7N+//6Lbr1ixAklJSfD19UV0dDR+/etfo7u7e5jSepaB7Puenh784Q9/QHx8PHQ6HdLT07Fhw4ZhTOs5srKysHDhQkREREChUGD9+vWXfExmZibGjx8PrVaLhIQErF69eshzeqqB7v+qqir86Ec/QmJiIpRKJR555JFhyemJBrrv161bhzlz5sBsNsNgMGDq1KnYuHHj8IT1MAPd9zt37sS0adMQHBwMX19fJCcn469//avL8rhNUfnwww+xZMkSPPnkk8jOzkZ6ejrmzZuH2traC27//vvv43e/+x2efPJJnDp1Cm+99RY+/PBD/P73vx/m5O5voPv+8ccfxxtvvIGXX34ZJ0+exAMPPIDbbrsNhw8fHubk7q+jowPp6el49dVXL2v7oqIi3HjjjZg1axZycnLwyCOP4P777+cb9iANdP9brVaYzWY8/vjjSE9PH+J0nm2g+z4rKwtz5szBV199hUOHDmHWrFlYuHAh33cGYaD73t/fHw8//DCysrJw6tQpPP7443j88cfx5ptvuiaQ5CYmT54sLV68uPdrh8MhRURESMuXL7/g9osXL5auu+66PvctWbJEmjZt2pDm9EQD3ffh4eHSK6+80ue+22+/XbrrrruGNKenAyB99tlnF93msccek1JTU/vcd+edd0rz5s0bwmTe4XL2/7fNmDFD+tWvfjVkebzJQPf9OSkpKdJTTz3l+kBeZLD7/rbbbpN+/OMfuySDWxxRsdlsOHToEGbPnt17n1KpxOzZs7Fnz54LPubqq6/GoUOHek9RFBYW4quvvsKCBQuGJbOnGMy+t1qt0Ol0fe7z9fXFzp07hzQrAXv27OnzswKAefPm9fuzIvJUTqcTbW1tCAoKEh3F6xw+fBi7d+/GjBkzXPJ8wmemvRz19fVwOBznzVYbFhaG06dPX/AxP/rRj1BfX49rrrkGkiTBbrfjgQce4KmfARrMvp83bx5eeOEFTJ8+HfHx8diyZQvWrVsHh8MxHJG9WnV19QV/Vq2trejq6oKvr6+gZETD6y9/+Qva29vxgx/8QHQUrxEVFYW6ujrY7XYsW7YM999/v0ue1y2OqAxGZmYmnnnmGbz22mvIzs7GunXr8OWXX+KPf/yj6Gge78UXX8SoUaOQnJwMjUaDhx9+GPfeey+USo/9dSMiGXn//ffx1FNP4aOPPkJoaKjoOF5jx44dOHjwIF5//XWsWLECH3zwgUue1y2OqISEhEClUqGmpqbP/TU1NbBYLBd8zP/+7//iJz/5SW+jS0tLQ0dHB37+85/jf/7nf/iheZkGs+/NZjPWr1+P7u5uNDQ0ICIiAr/73e8wcuTI4Yjs1SwWywV/VgaDgUdTyCusXbsW999/Pz7++OPzToPS0IqLiwNw9vO2pqYGy5Ytww9/+MMrfl63+LTWaDSYMGECtmzZ0nuf0+nEli1bMHXq1As+prOz87wyolKpAAASlze6bIPZ9+fodDpERkbCbrfj008/xS233DLUcb3e1KlT+/ysAGDTpk2X/FkReYIPPvgA9957Lz744APceOONouN4NafTCavV6pLncosjKgCwZMkSLFq0CBMnTsTkyZOxYsUKdHR04N577wUA3H333YiMjMTy5csBAAsXLsQLL7yAjIwMTJkyBfn5+fjf//1fLFy4sLew0OUZ6L7ft28fKioqMG7cOFRUVGDZsmVwOp147LHHRP413FJ7ezvy8/N7vy4qKkJOTg6CgoIQExODpUuXoqKiAmvWrAEAPPDAA3jllVfw2GOP4b777sPWrVvx0Ucf4csvvxT1V3BrA93/AJCTk9P72Lq6OuTk5ECj0SAlJWW447u1ge77999/H4sWLcKLL76IKVOmoLq6GsDZgfxGo1HI38FdDXTfv/rqq4iJiUFycjKAs5eK/+Uvf8Evf/lL1wRyybVDw+Tll1+WYmJiJI1GI02ePFnau3dv7/dmzJghLVq0qPfrnp4eadmyZVJ8fLyk0+mk6Oho6aGHHpKampqGP7gHGMi+z8zMlEaPHi1ptVopODhY+slPfiJVVFQISO3+tm3bJgE473Zufy9atEiaMWPGeY8ZN26cpNFopJEjR0qrVq0a9tyeYjD7/0LbjxgxYtizu7uB7vsZM2ZcdHu6fAPd9y+99JKUmpoq+fn5SQaDQcrIyJBee+01yeFwuCSPQpJ4HoSIiIjkyS3GqBAREZF3YlEhIiIi2WJRISIiItliUSEiIiLZYlEhIiIi2WJRISIiItliUSEiIiLZYlEhIiIi2WJRISIiItliUSEiIiLZYlEhIiIi2WJRISIiItn6/7VapZMRY1bSAAAAAElFTkSuQmCC",
      "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": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================== DoubleMLData 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",
      "Selection variable: s\n",
      "No. Observations: 8000\n",
      "\n",
      "------------------ DataFrame info    ------------------\n",
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 8000 entries, 0 to 7999\n",
      "Columns: 104 entries, X1 to s\n",
      "dtypes: float64(104)\n",
      "memory usage: 6.3 MB\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.DoubleMLData(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": 31,
   "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": 32,
   "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",
      "Selection variable: s\n",
      "No. Observations: 8000\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",
      "Learner ml_g_d0 RMSE: [[1.01990373]]\n",
      "Learner ml_g_d1 RMSE: [[1.19983954]]\n",
      "Learner ml_pi RMSE: [[0.37231324]]\n",
      "Learner ml_m RMSE: [[0.43374433]]\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": 33,
   "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": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABM20lEQVR4nO3deXxU5aE//s/sk2S2LDPZSEjIwr65gICKqHVBqdVb61UriNe2FryWcq2VXkVsq2hvpahVrP1aqL+K1gVt64IiCqjITgDZs5B9TyaTyUxmMjPn90dISiQsSWbmOTPzeb9e82oznMx8OETmw3Oe5zkKSZIkEBEREcmQUnQAIiIiojNhUSEiIiLZYlEhIiIi2WJRISIiItliUSEiIiLZYlEhIiIi2WJRISIiItliUSEiIiLZUosOMBSBQAA1NTUwGo1QKBSi4xAREdF5kCQJ7e3tyMjIgFJ59jGTiC4qNTU1yMrKEh2DiIiIBqGyshLDhg076zERXVSMRiOA7t+oyWQSnIaIiIjOh8PhQFZWVu/n+NlEdFHpudxjMplYVIiIiCLM+Uzb4GRaIiIiki2hRcXv9+PRRx9Fbm4u4uLikJeXh9/85jfgDZ2JiIgIEHzp5+mnn8aqVavw17/+FWPHjsWuXbswf/58mM1mPPDAAyKjERERkQwILSpbt27FTTfdhBtuuAEAkJOTg9dffx07duwQGYuIiGRGkiT4fD74/X7RUeg8qFQqqNXqoGwdIrSoTJ8+HS+//DKOHTuGwsJC7Nu3D19++SVWrFjR7/Eejwcej6f3a4fDEa6oREQkiNfrRW1tLVwul+goNADx8fFIT0+HVqsd0usILSoPP/wwHA4HRo0aBZVKBb/fjyeeeAJ33nlnv8cvX74cjz/+eJhTEhGRKIFAAGVlZVCpVMjIyIBWq+UGnzInSRK8Xi8aGxtRVlaGgoKCc27qdjZCi8qbb76J1157DWvXrsXYsWNRVFSERYsWISMjA/PmzTvt+CVLlmDx4sW9X/eswyYioujk9XoRCASQlZWF+Ph40XHoPMXFxUGj0aC8vBxerxd6vX7QryW0qPziF7/Aww8/jP/8z/8EAIwfPx7l5eVYvnx5v0VFp9NBp9OFOyYREQk2lH+RkxjB+jMT+ifvcrlO+42oVCoEAgFBiYiIiEhOhI6ozJkzB0888QSys7MxduxY7N27FytWrMA999wjMhYREclctd2N1g5v2N4vMUGLTEtc2N6P/k1oUXn++efx6KOPYsGCBWhoaEBGRgZ+8pOfYOnSpSJjERGRjFXb3bjqmU3o7Arf6Lteo8TG/7lCaFlZs2YNFi1aBLvdLiyDCEKLitFoxMqVK7Fy5UqRMYiIKIK0dnjR2RXAwln5YSkO1XY3Xvi8GK0d3gG93913342//vWvpz1/7bXXYv369Wf93pycHCxatAiLFi3qfe62227D7Nmzz/v9B0tuhSiib0pIRESxK9MSh9yUBNExzuq6667D6tWr+zw32EUhcXFxiIuLvctPnEZNREQUIjqdDmlpaX0eiYmJkCQJy5YtQ3Z2NnQ6HTIyMnpvHXPFFVegvLwcP//5z6FQKHr3jVmzZg0sFkvvay9btgyTJk3CX/7yF2RnZ8NgMGDBggXw+/343e9+h7S0NNhsNjzxxBN9Mq1YsQLjx49HQkICsrKysGDBAjidTgDApk2bMH/+fLS1tfW+97JlywB0b7r64IMPIjMzEwkJCZg6dSo2bdoU8nPIERWKWoGAhI++qcOmow3YW2HHdydl4L+vzOdmUUQk3DvvvIM//OEPeOONNzB27FjU1dVh3759AIB169Zh4sSJ+PGPf4wf/ehHZ32dkpISfPTRR1i/fj1KSkrw/e9/H6WlpSgsLMTmzZuxdetW3HPPPbj66qsxdepUAN3Lhp977jnk5uaitLQUCxYswEMPPYQXX3wR06dPx8qVK7F06VIcPXoUAGAwGAAA999/Pw4dOoQ33ngDGRkZePfdd3HdddfhwIEDKCgoCNm5YlGhqPX/vizFkx8eQXZSPDIseqzYcAzFDe343fcnQq9RiY5HRDHg/fff7/2g7/GrX/0Ker0eaWlpuPrqq6HRaJCdnY0pU6YAAJKSkqBSqWA0GpGWlnbW1w8EAvjLX/4Co9GIMWPGYNasWTh69Cg+/PBDKJVKjBw5Ek8//TQ+//zz3qJy6ryXnJwc/Pa3v8V9992HF198EVqtFmazGQqFos97V1RUYPXq1aioqEBGRgYA4MEHH8T69euxevVqPPnkk8E4Xf1iUaGodKjGgd+tP4obJ6TjzqnDAQCX5DbjxU0l8Pj24k93XSQ4IRHFglmzZmHVqlV9nktKSkJHRwdWrlyJESNG4LrrrsPs2bMxZ84cqNUD+1jOycmB0Wjs/To1NRUqlarPHmWpqaloaGjo/frTTz/F8uXLceTIETgcDvh8PnR2dsLlcp1x998DBw7A7/ejsLCwz/MejwfJyckDyjxQLCoUdTq7/PjZG3sxLDEOP7jo37dYmDoiGb6AhD9+XoydJ1pwcU6SwJREFAsSEhKQn59/2vNJSUk4evQoPv30U2zYsAELFizA//3f/2Hz5s3QaDTn/frfPlahUPT7XM9GqidOnMCNN96In/70p3jiiSeQlJSEL7/8Ev/1X/8Fr9d7xqLidDqhUqmwe/duqFR9R6S/PWIUbJxMS1HnlS/LcKK5AwuuyIdG1fdHfFpeMnJTEvDUR0cgSZKghERE3at45syZg+eeew6bNm3C119/jQMHDgAAtFot/H5/0N9z9+7dCAQCeOaZZ3DJJZegsLAQNTU1fY7p770nT54Mv9+PhoYG5Ofn93mc6/LUUHFEhaJKICBh7fYKTM9LQVbS6f8yUCoUuO2iLDy1/gg+O9KAq0anCkhJRMFQbXfL/n08Hg/q6ur6PKdWq/H+++/D7/dj6tSpiI+Px9/+9jfExcVh+PDuS9U5OTnYsmUL/vM//xM6nQ4pKSlD+j30yM/PR1dXF55//nnMmTMHX331FV566aU+x+Tk5MDpdGLjxo2YOHEi4uPjUVhYiDvvvBNz587FM888g8mTJ6OxsREbN27EhAkTcMMNNwQlX39YVCiqbC1pRrXdjR9dNuKMx0wYZsbYDBN+t/4orhxl4yogogiTmKCFXqPEC58Xh+099RolEhO0A/6+9evXIz09vc9zI0eOxFNPPYWnnnoKixcvht/vx/jx4/Gvf/2rd77Hr3/9a/zkJz9BXl4ePB5P0EaAJ06ciBUrVuDpp5/GkiVLcPnll2P58uWYO3du7zHTp0/Hfffdh9tuuw3Nzc147LHHsGzZMqxevRq//e1v8T//8z+orq5GSkoKLrnkEtx4441ByXYmCimCx78dDgfMZjPa2tpgMplExyEZ+O+1e7C7wo7ff3/CWQvIN9VteOLDw1i3YDouyE4MY0IiGojOzk6UlZUhNzcXer2+93ne60f+zvRnBwzs85sjKhQ17C4vPj5Yj1svGnbOUZIx6SYkJ2jx3t5qFhWiCJRpiWNxiBGcTEtR47291fBLEi4rsJ7zWKVSgWl5yfjnvhp0+cN3YzMiIhoYFhWKGu8V1WBylgXmuPNb2ndpfgrsri5sOdYY4mRERDRYLCoUFdpcXdhfZcekbMt5f8/w5ARkJ8Xj3b3VoQtGRERDwqJCUeHr0iYEJGBCpnlA3zcjPwUbDtWjvbMrRMmIKBgieN1HzArWnxmLCkWFL4ubkG7Ww2rUn/vgU8zIS4bHF8DnR3n5h0iOenZZdblcgpPQQPX8mQ1kp93+cNUPRYUvjjVhbMbARlMAINmgQ3ZSPDYfbcR3J2aEIBkRDYVKpYLFYum9V018fDz3PpI5SZLgcrnQ0NAAi8Vy2pb7A8WiQhGvssWF8hYX/uOCYYP6/gnDzNh0rAGBgASlkn8BEslNzxbtp95Yj+TPYrEEZXt9FhWKeF8WN0GpAMZkDG7Tv0lZFry/vxaHah0YN8A5LkQUegqFAunp6bDZbOjq4nyySKDRaIY8ktKDRYUi3hfHG5FnMyBBN7gf55GpRug1Smw+1siiQiRjKpUqaB9+FDk4mZYiWiAg4aviZowfxPyUHmqVEmMzzNh0lMPKRERyw6JCEa2suQNt7i6MTDMO6XUmDrNgT7kdDi5TJiKSFRYVimj7Ku0AgBFWw5BeZ+IwM/yShK3FTUFIRUREwcKiQhFtf1UbMix6GAY5P6WHzaRHhkWPL46zqBARyQmLCkW0vZWtyE0Z2mhKj1FpJmwvawnKaxERUXCwqFDE8voCOFzTjnxrQlBeb1SaEcUNTjQ7PUF5PSIiGjoWFYpYR+va4fUHkDfE+Sk9xqR378Oy80RrUF6PiIiGjkWFIlZRlR0qpQLDk4MzopJs0MFm1GEHL/8QEckGiwpFrP2VdgxPiodWHbwf45FpRmwrbQ7a6xER0dCwqFDEKqq0IzclOKMpPUanmXC41sH9VIiIZIJFhSKS0+NDcYMTebbgzE/pMSrdCAnAbs5TISKSBRYVikgHq9sgAUGbSNsjzaRHUryGy5SJiGSCRYUi0pG6dqhVCmRY9EF9XYVCgZFpJmwv4zwVIiI5YFGhiHSkrh3DLHFQK4P/I1yYasA31W3w+PxBf20iIhoYoUUlJycHCoXitMfChQtFxqIIcKTOgczE+JC8dkGqEV1+CQdrHCF5fSIiOn9Ci8rOnTtRW1vb+9iwYQMA4NZbbxUZi2ROkiQcq2tHdmJcSF5/eFI8tCol9lbYQ/L6RER0/oQWFavVirS0tN7H+++/j7y8PMycOVNkLJK5arsbHV4/hiWFZkRFrVIi15qAPRVc+UNEJNrQbjkbRF6vF3/729+wePFiKBSKfo/xeDzweP59HxaHg0PzsehYfTsAIDtERQUACmwG7C5nUSEiEk02k2nfe+892O123H333Wc8Zvny5TCbzb2PrKys8AUk2ThS1454rQrJCdqQvUe+zYDatk7UOzpD9h5ERHRusikqr7zyCq6//npkZGSc8ZglS5agra2t91FZWRnGhCQXx+rakZUYf8aRt2AosBkBAHs4qkJEJJQsLv2Ul5fj008/xbp16856nE6ng06nC1MqkqvDde0YFqKJtD2SErSwGnTYW2nH9ePTQ/peRER0ZrIYUVm9ejVsNhtuuOEG0VFI5rr8AZQ0OJEVwvkpPfI5T4WISDjhRSUQCGD16tWYN28e1GpZDPCQjJ1o6oAvIIWtqByoaoPXFwj5exERUf+EF5VPP/0UFRUVuOeee0RHoQhwpK57xU9WiC/9AN0rf7z+AI6efE8iIgo/4UMY11xzDSRJEh2DIsSx+nYkxmtg1GtC/l7DkxOgUipQVGXH+GHmkL8fERGdTviICtFAlDQ6kWEJ/WgKAGjVSuQkx2NfpT0s70dERKdjUaGIUtzgRLo5uHdMPpvcFAOKuJU+EZEwLCoUMQIBCeXNLqSbwzOiAgD5tgSUNDrR3tkVtvckIqJ/Y1GhiFFtd8PjCyDDEr4RlTyrARKAA9VtYXtPIiL6NxYVihhlTR0AENYRlQxzHOI0KuyrZFEhIhKBRYUiRmmjExqVAlZD+HYnVioVGGFNwL5KbvxGRCQCiwpFjNKmDqSZ9VAqQ3ePn/7kWQ0o4ogKEZEQLCoUMUoanUgzhW9+So88qwF1Dt5JmYhIBBYVihgljR1hnZ/SI8+aAAAo4n4qRERhx6JCEcHl9aGurTOsK356JCVokRivwf4qe9jfm4go1rGoUEQQseKnh0KhwAirgSMqREQCsKhQROgpKhkCigoAjEhJwP6qNgQCvC8VEVE4sahQRCht7IBJr4ZBL+Y+mvk2A9o7fTjR3CHk/YmIYhWLCkWE0kYn0sN0M8L+jEgxAAD2V3GZMhFROLGoUEQoaewQsjS5h0GvRrpZz3kqRERhxqJCEaG8uXuzN5FGpCRgH1f+EBGFFYsKyZ7d5YWj04dUo+CiYjXgUI0DXf6A0BxERLGERYVkr7zZBQBINYXvHj/9ybcZ4PEFcLSuXWgOIqJYwqJCslfe0lNUxI6o5CQnQKkAL/8QEYURiwrJXkVzB4x6NRJ0YpYm99CqlRieHI99nFBLRBQ2LCoke+XNLuGjKT1GpBiwt8IuOgYRUcxgUSHZO9HcAZtR7PyUHnk2A4obnHB6fKKjEBHFBBYVkj05jajkWw2QABzgxm9ERGHBokKy1tnlR0O7R/iKnx6ZljjEaVTc+I2IKExYVEjWKnpW/AjeQ6WHUqnACGsCiipbRUchIooJLCokaydO3jU5VfCutKfKsxo4okJEFCYsKiRrFS0u6NRKWOI0oqP0yrcaUO/woK6tU3QUIqKox6JCstYzkVahUIiO0ivP1n0nZV7+ISIKPRYVkrVyGS1N7pGUoEVyghZFlVz5Q0QUaiwqJGsnZLQ0+VR5VgP2VnBEhYgo1FhUSLZ8/gBq7G7ZLE0+VZ7NgP1VbfDxTspERCHFokKyVdvWCV9Agk0mS5NPVZhqgLvLjyO8kzIRUUixqJBsVZ7cQ0Vuc1SA7nv+qJQKXv4hIgoxFhWSrapWNxQAUmRYVLRqJXKT47GHNygkIgop4UWluroaP/zhD5GcnIy4uDiMHz8eu3btEh2LZKCy1YWkBC00KuE/pv3Ktxmxu5wjKkREoST0E6C1tRUzZsyARqPBRx99hEOHDuGZZ55BYmKiyFgkE5UtLlmOpvQoTDWgosWFJqdHdBQioqilFvnmTz/9NLKysrB69ere53JzcwUmIjmpaHHBapBvUSlINQIA9pS34pqxaYLTEBFFJ6EjKv/85z9x0UUX4dZbb4XNZsPkyZPx5z//WWQkkpHKVrcsJ9L2SE7QIilBy3kqREQhJLSolJaWYtWqVSgoKMDHH3+Mn/70p3jggQfw17/+td/jPR4PHA5HnwdFp84uPxrbPbDKuKgoFArk2wzYXd4iOgoRUdQSWlQCgQAuuOACPPnkk5g8eTJ+/OMf40c/+hFeeumlfo9fvnw5zGZz7yMrKyvMiSlcqu1uAPJcmnyqApsBB6ra0MWN34iIQkJoUUlPT8eYMWP6PDd69GhUVFT0e/ySJUvQ1tbW+6isrAxHTBKgZw8VOY+oAMDIVCM6fQF8U837/hARhYLQybQzZszA0aNH+zx37NgxDB8+vN/jdToddDp5f3BRcFS2uqFUAEkJ8v7zzk1JgE6txM4TLZiczdVqRETBJnRE5ec//zm2bduGJ598EsXFxVi7di1efvllLFy4UGQskoGqVhesRh1USoXoKGelVimRbzNgRxnnqRARhYLQonLxxRfj3Xffxeuvv45x48bhN7/5DVauXIk777xTZCySgaoWt6yXJp9qZJoRO0+0IhCQREchIoo6Qi/9AMCNN96IG2+8UXQMkpnyFpfs56f0GJ1mwro91ShudKLw5N4qREQUHPLcm5xiXlWLCykRMqKSb+u+QeF2Xv4hIgo6FhWSHafHB7u7CzaTXnSU86LXqJCbkoCdLCpEREHHokKyU9XavTRZ7nuonGpUmhE7ylogSZynQkQUTCwqJDuVLd2bvUXKpR+ge0JtnaMTVa1u0VGIiKIKiwrJTlWrCxqVApZ4jego521UqgkAOE+FiCjIWFRIdqpb3Ugx6KBUyHsPlVMZ9GrkJMfj65Jm0VGIiKIKiwrJTrXdHVGXfXqMyTDjq+ImzlMhIgoiFhWSncpWN1IMWtExBmxchgl1jk6UNXWIjkJEFDVYVEh2qlsjZw+VU41KM0GlVGArL/8QEQUNiwrJisvrQ6urK2J2pT1VnFaFfKsBXxU3iY5CRBQ1WFRIVmrskbc0+VRjMkz4uqSZ9/0hIgoSFhWSlcqT+5BE4ogK0D1Pxe7uwuE6h+goRERRgUWFZKW61Q2lAkiMj7zJtABQkGqEVqXE1mLOUyEiCgYWFZKVnqXJKmXk7KFyKo1KiVFpRnxxvFF0FCKiqMCiQrLSs9lbJJswzIJtpS1we/2ioxARRTwWFZKVylYXkiNwD5VTTcqywOsPYFspL/8QEQ0ViwrJSnWrG9YIH1HJsOhhM+qw6WiD6ChERBGPRYVkw+Pzo7HdE/GXfhQKBSYMs+Dzo5ynQkQ0VCwqJBu19k5IiNylyaealGVBRYuL2+kTEQ0RiwrJRnWEb/Z2qrEZJmhUCl7+ISIaIhYVko3qk5u9RfpkWgDQa1QYlWbC50dYVIiIhoJFhWSjyu5GUoIWGlV0/FhOyupeptzh8YmOQkQUsaLjE4GiQvceKpE/mtLjouGJ8PoD2HyMk2qJiAaLRYVko6rVheSEyJ+f0sNm0iMnOR4fH6wTHYWIKGKxqJBsdG+fHz0jKgBw4fBEfHa4AV3+gOgoREQRiUWFZMEfkFDX1hkVK35OdVFOEto9PmwvbREdhYgoIrGokCw0tnvgC0hRV1SGJ8XDatTx8g8R0SCxqJAs9O6hEgWbvZ1KoVDgwuGJ+ORQHQIBSXQcIqKIw6JCsvDvzd6ia44KAFw8PBH1Dg+KquyioxARRRwWFZKFGrsbCVoV4rVq0VGCblSaCYnxGnywv1Z0FCKiiMOiQrJQ3epGcpTNT+mhVCpwcU4S3t9fw8s/REQDxKJCshCNS5NPNS0vGfUOD3aVt4qOQkQUUVhUSBaqWl1RO6ICAIWpRiQnaPH+/hrRUYiIIgqLCslCrT369lA5lVKhwNTcJHxwoBZ+Xv4hIjpvQovKsmXLoFAo+jxGjRolMhIJ4OjsQrvHF9WXfoDuyz/NTi+2lTaLjkJEFDGEL7EYO3YsPv30096v1WrhkSjManqXJkfviAoA5FkNSDXp8I+iaszITxEdh4goIgi/9KNWq5GWltb7SEnhX+Cxpro1NoqKQqHAjLwUfHigDp1dftFxiIgigvCicvz4cWRkZGDEiBG48847UVFRIToShVmN3Q21UgFLvEZ0lJC7tCAFTo8Pnx6uFx2FiCgiCC0qU6dOxZo1a7B+/XqsWrUKZWVluOyyy9De3t7v8R6PBw6Ho8+DIl+V3Y1kgxZKhUJ0lJBLN8ehwGbAuj3VoqMQEUUEoUXl+uuvx6233ooJEybg2muvxYcffgi73Y4333yz3+OXL18Os9nc+8jKygpzYgqFmihf8fNtM/JTsPlYI5qdHtFRiIhkT/iln1NZLBYUFhaiuLi4319fsmQJ2traeh+VlZVhTkihUNXqQlJCdK/4OdW0vGQAwL/2cU8VIqJzkVVRcTqdKCkpQXp6er+/rtPpYDKZ+jwo8tXY3TE1omLSazBpmAXv8PIPEdE5CS0qDz74IDZv3owTJ05g69atuPnmm6FSqXD77beLjEVh5PUF0ODwxFRRAYDLC604UN2GY/X9z8ciIqJuQotKVVUVbr/9dowcORI/+MEPkJycjG3btsFqtYqMRWFU7+iEBET9Zm/fdkG2BUa9Gu/srhIdhYhI1oTurvbGG2+IfHuSgeoY2ezt29QqJabnpeCdPVX4xbUjoVbJ6iosEZFs8G9HEqpns7fkGBtRAYCZhVY0Ob344niT6ChERLLFokJC1djdMMdpoFOrREcJu5zkeGQnxeOt3Vy9RkR0JiwqJFRNmzvm5qf0UCgUuLzAig2H6mF3eUXHISKSJRYVEqqq1R1Te6h824z8ZPgDEv61v1Z0FCIiWWJRIaGqW2NrD5Vvs8RrMSnLgrd28fIPEVF/WFRIGEmSTl76id2iAnTvqbK/qg3HuacKEdFpWFRImFZXFzq7AjG54udUF2QnwqBT4+093FOFiOjbWFRImJoY3UPl2zQqJabnJePdPdXwByTRcYiIZIVFhYSpamVR6XF5oRUN7R58cbxRdBQiIllhUSFhauxuaFVKmPRCN0iWhREpCci0xOHdvbxRIRHRqVhUSJgauxspRi0UCoXoKMIpFApcmp+Cjw/WocPjEx2HiEg2BlVUSktLg52DYlC13Y3kBF726TEjPwWdXQGs/6ZOdBQiItkYVFHJz8/HrFmz8Le//Q2dnZ3BzkQxoqo1dnel7Y/VqMPodCMv/xARnWJQRWXPnj2YMGECFi9ejLS0NPzkJz/Bjh07gp2NolxNmxvJnEjbx6X5VmwtaUK9g/8AICICBllUJk2ahGeffRY1NTX4y1/+gtraWlx66aUYN24cVqxYgcZGrlygs+vs8qPZ6eWIyrdMzU2CSqnAP4tqREchIpKFIU2mVavVuOWWW/DWW2/h6aefRnFxMR588EFkZWVh7ty5qK3l/Uuof7Vt3SMGnKPSV4JOjclZifjnPhYVIiJgiEVl165dWLBgAdLT07FixQo8+OCDKCkpwYYNG1BTU4ObbropWDkpylSf3EPFamRR+bZLRiTjQHUbyps7REchIhJuUBtYrFixAqtXr8bRo0cxe/ZsvPrqq5g9ezaUyu7ek5ubizVr1iAnJyeYWSmKVNtdUABIjuE7J5/J5GwL9Bol3t9fi4Wz8kXHISISalAjKqtWrcIdd9yB8vJyvPfee7jxxht7S0oPm82GV155JSghKfpU2zuRmKCFWsWtfL5Nr1Hhgmxe/iEiAgY5orJhwwZkZ2efVk4kSUJlZSWys7Oh1Woxb968oISk6FPNpclndcmIZKzYcAzFDe3ItxlFxyEiEmZQ/5zNy8tDU1PTac+3tLQgNzd3yKEo+lW1ujiR9iwmDrMgXqvCv/ZxQjoRxbZBFRVJ6v8Or06nE3q9fkiBKDZU2zmicjZatRIXZifigwMsKkQU2wZ06Wfx4sUAuu9LsnTpUsTHx/f+mt/vx/bt2zFp0qSgBqTo4w9IqGvrRMoYjqiczcW5SVix4RhKGp3IsxpExyEiEmJARWXv3r0AukdUDhw4AK323/8i1mq1mDhxIh588MHgJqSo09jugS8gIYWXfs5qwjAzdGolPj5YhwVXcPUPEcWmARWVzz//HAAwf/58PPvsszCZTCEJRdGt2u4CAKRwD5Wz0qlVmDjMgvXfsKgQUewa1ByV1atXs6TQoFWd3OyNc1TO7aKcROyvakNdG+/9Q0Sx6bxHVG655RasWbMGJpMJt9xyy1mPXbdu3ZCDUfSqsXciQadCvHZQq+NjyuTsRKiUCnxyqA5zp+WIjkNEFHbn/UlhNpuhUCh6/z/RYFXbXUjhXZPPi0Gnxth0E9Z/w6JCRLHpvIvK6tWr+/3/RANV3ermRNoBuCgnEX/dWo42VxfM8RrRcYiIwmpQc1TcbjdcLlfv1+Xl5Vi5ciU++eSToAWj6FXV6kYy56ectwuyE+GXJGw53ig6ChFR2A2qqNx000149dVXAQB2ux1TpkzBM888g5tuugmrVq0KakCKLpIkndzsjSMq5yvZoMPw5Hh8dqRBdBQiorAbVFHZs2cPLrvsMgDA22+/jbS0NJSXl+PVV1/Fc889F9SAFF0cbh9cXj+LygBNyrLg86MN8Af63xWaiChaDaqouFwuGI3dN0r75JNPcMstt0CpVOKSSy5BeXl5UANSdKk6uYeK1chLPwMxOSsRdlcXiirtoqMQEYXVoIpKfn4+3nvvPVRWVuLjjz/GNddcAwBoaGjg/ip0VjX27v1AOKIyMAU2A4x6NT7n5R8iijGDKipLly7Fgw8+iJycHEydOhXTpk0D0D26Mnny5KAGpOhS3eqCRqWAKY6rVwZCqVRgQqYZG4/Ui45CRBRWgyoq3//+91FRUYFdu3Zh/fr1vc9fddVV+MMf/jCoIE899RQUCgUWLVo0qO+nyFDV2j2RVnlyTx46f5OzE3G4tp271BJRTBn01qBpaWlIS0vr89yUKVMG9Vo7d+7En/70J0yYMGGwcShCVLW6YeVln0GZMMwMBYDNxxpw28XZouMQEYXFoEZUOjo68Oijj2L69OnIz8/HiBEj+jwGwul04s4778Sf//xnJCYmDiYORZCqVhdvRjhIRr0GedYEfHG8SXQUIqKwGdSIyr333ovNmzfjrrvuQnp6eu/W+oOxcOFC3HDDDbj66qvx29/+9qzHejweeDye3q8dDseg35fEqGp1Y0wGb8EwWOMyzdh8rBGBgASlkpfPiCj6DaqofPTRR/jggw8wY8aMIb35G2+8gT179mDnzp3ndfzy5cvx+OOPD+k9SZwOjw92dxesHFEZtPHDLHivqAaHah0Yl8nCR0TRb1CXfhITE5GUlDSkN66srMTPfvYzvPbaa9Dr9ef1PUuWLEFbW1vvo7KyckgZKLyq7W4A4ByVISi0GaDXKHn5h4hixqCKym9+8xssXbq0z/1+Bmr37t1oaGjABRdcALVaDbVajc2bN+O5556DWq2G3+8/7Xt0Oh1MJlOfB0WOqtaezd5YVAZLrVJidJoJW47xvj9EFBsGdennmWeeQUlJCVJTU5GTkwONpu+eGHv27Dnna1x11VU4cOBAn+fmz5+PUaNG4Ze//CVUKtVgopGMVbW6oVYqYOEdgIdk/DAzXt9RAbfXjzgt/zshoug2qKLyve99b8hvbDQaMW7cuD7PJSQkIDk5+bTnKTpUcw+VoBifacarfgk7TrRgZqFVdBwiopAaVFF57LHHgp2DYkBVqxspvMfPkGVa4pCcoMWXxxtZVIgo6g16wze73Y63334bJSUl+MUvfoGkpCTs2bMHqampyMzMHNRrbtq0abBxKAJUtLiQksD5KUOlUCgwOt2ErSXNoqMQEYXcoCbT7t+/H4WFhXj66afx+9//Hna7HQCwbt06LFmyJJj5KIpU292cSBskYzJMOFTjQJurS3QUIqKQGlRRWbx4Me6++24cP368z9Li2bNnY8uWLUELR9HD7fWjpcPLohIkY9NNkABsL+OoChFFt0EVlZ07d+InP/nJac9nZmairq5uyKEo+lTbu5cmp3APlaCwmfSwGnXYVtoiOgoRUUgNqqjodLp+t68/duwYrFZO7qPTVbWe3OyNIypBMybdhK0l3PiNiKLboIrKd7/7Xfz6179GV1f39XGFQoGKigr88pe/xH/8x38ENSBFh6pWN1RKBRLjueonWEanm3Ckrh2tHV7RUYiIQmZQReWZZ56B0+mE1WqF2+3GzJkzkZ+fD6PRiCeeeCLYGSkKVLW6kWLQQsUb6QXN2IzunZm3l/HyDxFFr0EtTzabzdiwYQO++uor7Nu3D06nExdccAGuvvrqYOejKFFtd3N+SpClGHRINemwrbQZ141LEx2HiCgkBlxUAoEA1qxZg3Xr1uHEiRNQKBTIzc1FWloaJEmCgruOUj/KmztYVEKA81SIKNoN6NKPJEn47ne/i3vvvRfV1dUYP348xo4di/Lyctx99924+eabQ5WTIlxViws2TqQNulFpJhyrd8Lu4jwVIopOAxpRWbNmDbZs2YKNGzdi1qxZfX7ts88+w/e+9z28+uqrmDt3blBDUmTr8PjQ4uriip8QGJVmBADsOtGKq8ekCk5DRBR8AxpRef311/GrX/3qtJICAFdeeSUefvhhvPbaa0ELR9GhZ2myzag/x5E0UFajDskJWuw8wQm1RBSdBlRU9u/fj+uuu+6Mv3799ddj3759Qw5F0aWypXuzN46oBJ9CoUBhmhHbuEMtEUWpARWVlpYWpKaeeXg5NTUVra2tQw5F0aWy1QWNSgFLvEZ0lKg0Os2Ig9UOuLw+0VGIiIJuQEXF7/dDrT7ztBaVSgWfj39ZUl+VLW7YjHoouSIsJEalmeALSCiqsIuOQkQUdAOaTCtJEu6++27odP0P4Xs8nqCEouhS0dLByz4hlJkYB4NOjR0nWjA9P0V0HCKioBpQUZk3b945j+GKH/q2ihYXhicniI4RtZQKBUamGrGDO9QSURQaUFFZvXp1qHJQlJIkCVWtblyckyQ6SlQbmWbEur1V6PIHoFEN6s4YRESyxL/RKKRaXV1wef289BNio9KM6OwK4FDN6Xc1JyKKZCwqFFI9S5O5h0po5aYkQKNSYFc5V90RUXRhUaGQqmzlHirhoFYpkWc1YBc3fiOiKMOiQiFV2eJGgk4Fg25QN+qmAShMNWJ3eSskSRIdhYgoaFhUKKQqWly87BMmhalGNLR7em9ZQEQUDVhUKKQqW1ywGnjZJxwKUg0AgD0VnKdCRNGDRYVCqqLFxfkpYWLSa5BpicOuEywqRBQ9WFQoZPwBCTV2N2wmFpVwybcZeCdlIooqLCoUMjV2N3wBCamcoxI2I1ONOFbfjvbOLtFRiIiCgkWFQqbi5B4qqSYWlXApTDUiIAFFlXbRUYiIgoJFhULmRHMHlAogxagVHSVmZFj0MOjU2M2N34goSrCoUMhUNLtgNeqhVvLHLFwUCgUKbAbsYVEhoijBTxAKmfJmF1K54ifs8m0G7K2wIxDgxm9EFPlYVChkypo7uOJHgMJUI9o9PpQ0OkVHISIaMhYVCglJklDR7OJEWgHyrAYoFdz4jYiiA4sKhUST0wt3l59FRYA4rQpZSfGcUEtEUYFFhUKioqUDAJcmi1JgM7CoEFFUEFpUVq1ahQkTJsBkMsFkMmHatGn46KOPREaiIDnR1L2Hio2TaYUosBlR0tiBNhc3fiOiyCa0qAwbNgxPPfUUdu/ejV27duHKK6/ETTfdhIMHD4qMRUFQ3uJCYrwGeo1KdJSY1HODwr2VHFUhosgmtKjMmTMHs2fPRkFBAQoLC/HEE0/AYDBg27ZtImNREFQ0d/Cyj0BpJj1MejX2VNhFRyEiGhK16AA9/H4/3nrrLXR0dGDatGn9HuPxeODxeHq/djgc4YpHA3Si2cXLPgIpFArkc+M3IooCwifTHjhwAAaDATqdDvfddx/effddjBkzpt9jly9fDrPZ3PvIysoKc1o6X+UcURGuwGZEUSU3fiOiyCa8qIwcORJFRUXYvn07fvrTn2LevHk4dOhQv8cuWbIEbW1tvY/Kysowp6Xz0d7ZhVZXF4uKYAWpBjg9Phxv4MZvRBS5hF/60Wq1yM/PBwBceOGF2LlzJ5599ln86U9/Ou1YnU4HnY6XE+SuvLnnrsn8sxKpZ+O3vRWtGJlmFB2HiGhQhI+ofFsgEOgzD4UiT1lT9x4qaaY4wUlim16jQnZSPHeoJaKIJnREZcmSJbj++uuRnZ2N9vZ2rF27Fps2bcLHH38sMhYNUVlTB0x6NQx64QN2MS+fG78RUYQT+knS0NCAuXPnora2FmazGRMmTMDHH3+M73znOyJj0RCVNXUg3cLRFDkosBnx6eEGtLm7YI7TiI5DRDRgQovKK6+8IvLtKURKGp1I40RaWSiwdW/8VlRpx8xCq+A0REQDJ7s5KhTZJEnqHlExs6jIQZpZD6Nezf1UiChisahQUDV3eNHe6UO6mZd+5KB34zdOqCWiCMWiQkHVs+KHIyrykW81oKiCG78RUWRiUaGgKmvsgALgZm8yUphqRLvHh+JGbvxGRJGHRYWCqrSpA1ajDlo1f7TkIt/WvfEb56kQUSTipwkFVVkTV/zIDTd+I6JIxqJCQVXa2IE0zk+RHW78RkSRikWFgsYfkFDe7OJEWhkqsBlR0tgBu8srOgoR0YCwqFDQ1Njd8PoDSOPSZNkpTO2+KeHeSrvYIEREA8SiQkHDpcnylWrSwaRXYy8v/xBRhGFRoaApbXRCrVTAatCJjkLfolAoUGAzYjcn1BJRhGFRoaApbnQiwxIHpVIhOgr1Iz+1e+M3Pzd+I6IIwqJCQXO83snLPjJWmGpEh9ePY/XtoqMQEZ03FhUKmuMNTgxL5ERaucqzJkClVHCZMhFFFBYVCoqWDi9aOrzItMSLjkJnoFOrkJscz6JCRBGFRYWCorih+z4yHFGRt/xUI3aeaBEdg4jovLGoUFAUNzihVIC70srcyFQjqlrdaHB0io5CRHReWFQoKI43tCPNrIdGxR8pOevZ+I2Xf4goUvBThYLieL0TmRZe9pG7pAQtrEYddrGoEFGEYFGhoDje0M6iEiEKbQbsKuc8FSKKDCwqNGTtnV2od3iQmcgVP5GgMM2Ig9UOdHb5RUchIjonFhUasp4VPxxRiQyFqUb4AhL28QaFRBQBWFRoyI43OKEAkGHhip9IkJ0YjziNivNUiCgisKjQkJU0OGEz6aBTq0RHofOgVCpQmGrAjjLOUyEi+WNRoSE7Vt+ODDMv+0SSkWkm7C5v5Q0KiUj2WFRoyA7XtiMriRNpI8moNCOcHh+O1DlERyEiOisWFRoSu8uLOkcnhiezqESSPKsBGpUCO3n5h4hkjkWFhuRwbTsAIJsjKhFFq1ZihNWAHbzvDxHJHIsKDcmROgc0KgXSOUcl4oxMNWJHWQskifNUiEi+WFRoSA7XOpCVGA+VUiE6Cg3QyDQjmpxelDe7REchIjojFhUakkM1Dk6kjVAjU41QALz8Q0SyxqJCg+bzB3C8wcn5KREqQadGdnI8J9QSkayxqNCgnWjugMcX4IqfCDYqzYRtpc2iYxARnRGLCg0aV/xEvjHpJlS2ulFtd4uOQkTUL6FFZfny5bj44othNBphs9nwve99D0ePHhUZiQbgcK0DyQlaGPUa0VFokEanGwEA20o4qkJE8iS0qGzevBkLFy7Etm3bsGHDBnR1deGaa65BR0eHyFh0no7UOjiaEuGMeg2GJ8fz8g8RyZZa5JuvX7++z9dr1qyBzWbD7t27cfnllwtKRefrUG07puQmiY5BQzQ6zYSvWVSISKZkNUelra0NAJCU1P+Hn8fjgcPh6PMgMZqdHtQ5OpHDibQRb0yGCVWtblS1cj8VIpIf2RSVQCCARYsWYcaMGRg3bly/xyxfvhxms7n3kZWVFeaU1GN/dXepHGE1CE5CQzU6zQQFgG2lXKZMRPIjm6KycOFCfPPNN3jjjTfOeMySJUvQ1tbW+6isrAxjQjrVgao2GPVq2Iw60VFoiAx6NeepEJFsCZ2j0uP+++/H+++/jy1btmDYsGFnPE6n00Gn4wejHOyrsiM3JQEKBbfOjwaj0k3YWtIkOgYR0WmEjqhIkoT7778f7777Lj777DPk5uaKjEPnSZIk7Ku0Y0RKgugoFCRjM0yosXeivJkr7ohIXoSOqCxcuBBr167FP/7xDxiNRtTV1QEAzGYz4uJ4N165qnd40OT0cn5KFBmTboJSAXxZ3IThySygRCQfQkdUVq1ahba2NlxxxRVIT0/vffz9738XGYvOYV+VHQA4ohJF4rVqFNiM+OI4L/8QkbwIHVGRJEnk29MgHahqQ2K8BkkJWtFRKIjGZZqw4VA9/AEJKiXnHhGRPMhm1Q9FDk6kjU7jMs1wdPrwzcml50REcsCiQgMiSRL2V7VxfkoUyrcZEKdR4ctiXv4hIvlgUaEBqWxxo83dhTwr56dEG7VSiTHpJs5TISJZYVGhAdlb2QoAGJHCEZVoNC7ThN3lLXB7/aKjEBEBYFGhAdp5ogWZljiY4jSio1AIjM+0oMsvYVsZd6klInlgUaEB2VHWgsJUo+gYFCIZFj2sRh02H20UHYWICACLCg1Am6sLx+udGJnGohKtFAoFJg6z4LMjDaKjEBEBYFGhAdhT0QoJwCgWlag2OcuCihYXypq4nT4RiceiQudt54kWJMZreMfkKDcmwwS1SoHPOapCRDLAokLnbeeJ7vkp3Ogtuuk1KoxJN+HzoywqRCQeiwqdF4/Pj32VbZyfEiMmZVmwvbQFLq9PdBQiinEsKnReDlS1wesPYCRX/MSESVkWeP0BfF3CZcpEJBaLCp2XHSdaoNcoMTyZO9LGgnRzHNLNenx6mJd/iEgsFhU6L18VN2FUmpF31Y0hF2QnYsOhOgQCvMs5EYnDokLn1Nnlx86yVozPtIiOQmF0UU4impxeFFXZRUchohjGokLntKOsBV5/AOMzzaKjUBgV2owwx2nwycF60VGIKIaxqNA5fVnchKQELYYlxomOQmGkVCpwQbYFHx+sEx2FiGIYiwqd0+ZjjRiXYeL+KTHoouFJKGvqQHGDU3QUIopRLCp0Vo3tHhyta8f4YRbRUUiAcZlm6DVKfHKIoypEJAaLCp3VV8VNAIBxGSbBSUgErVqJCcMs+OgAiwoRicGiQme15XgjcpLjYYnXio5CglySm4QD1W2oaHaJjkJEMYhFhc7IH5Cw6WgjV/vEuMnZidCplXj/QI3oKEQUg1hU6Ix2l7eipcOLi3OSREchgfQaFSZnW/CvfSwqRBR+LCp0Ruu/qUNivAZ5NoPoKCTYtBEpOFzbjtJGrv4hovBiUaF+SZKE9QdrceHwJCi5LDnmTcqyIE6jwvv7a0VHIaIYw6JC/TpY40CNvRMX5ySKjkIyoFUrceHwRPyTl3+IKMxYVKhfnxysQ4JWhTHpXJZM3ablJaO4wYmDNW2ioxBRDGFRoX6tP1iHSdmJUKv4I0LdJgwzwxynwTu7q0VHIaIYwk8hOk1xQzuO1Tsxhat96BRqpRIz8pLxXlE1uvwB0XGIKEawqNBp3tlTjQRd95JUolNdXmhFS4cXm482io5CRDGCRYX68AckrNtThWkjkqHhZR/6luHJCchJjsfbe6pERyGiGMFPIupja0kT6h0eXF5gFR2FZOqyAis2Hq6H3eUVHYWIYgCLCvXxzu4qZJj1yOcmb3QGM/JTEJCAdXs4qZaIQo9FhXo5PT6sP1iHSwusUHCTNzoDc5wGF+ck4rXt5ZAkSXQcIopyQovKli1bMGfOHGRkZEChUOC9994TGSfm/WtfDTxdAVxWkCI6Csnc1aNTUdLYgR1lLaKjEFGUE1pUOjo6MHHiRLzwwgsiYxC6t8xf/VUZLhyeiBSDTnQckrkx6SZkmPV4bXu56ChEFOXUIt/8+uuvx/XXXy8yAp30dWkzjtU78b+zR4uOQhFAoVDgqtGpeH1HBZbO8bDcElHIRNQcFY/HA4fD0edBwbH6yxPISozD2AxumU/n5/ICKxQK4M1dlaKjEFEUi6iisnz5cpjN5t5HVlaW6EhRobLFhU8P1+PasWmcREvnzaBXY3peCl7dWs6daokoZCKqqCxZsgRtbW29j8pK/ksuGFZ/dQIJOjVm5HMSLQ3M7PHpqHN04sMDtaKjEFGUiqiiotPpYDKZ+jxoaBrbPVi7vRzXjEmFXqMSHYciTHZSPCYMM+PlLaVcqkxEIRFRRYWC7+UtJVAqFbh+fLroKBShrh+XjoM1Di5VJqKQELrqx+l0ori4uPfrsrIyFBUVISkpCdnZ2QKTxYbGdg/+v6/LMXt8Ogw6oT8KFMEmDjNjWGIcXt5SiqkjkkXHIaIoI3REZdeuXZg8eTImT54MAFi8eDEmT56MpUuXiowVM3pHU8ZxNIUGT6FQ4Ibx6dh4pAGHa7kSj4iCS2hRueKKKyBJ0mmPNWvWiIwVE6rtbrz6dTmuG5cGg56jKTQ0lxakwGbU4fnPjouOQkRRhnNUYtTTHx1GvFaFORMyREehKKBWKvHdiRn46EAdjte3i45DRFGERSUG7S5vxT/31eLWi7K40oeC5vJCK5INWvzxs+JzH0xEdJ5YVGJMICDh1+8fRE5yPGYWWEXHoSiiUSkxZ2IG/rW/BsUNHFUhouBgUYkxb+6qxL7KNtx1yXAoldyFloJr1kgbUgw6PP3RUdFRiChKsKjEkIb2Tjzx4WHMLLRiTIZZdByKQhqVErdelIUNh+ux6wT3VSGioWNRiSGP//MQlAoF7pzKPWoodKbnJSMnOR5PfniYu9US0ZCxqMSIDYfq8cGBWtx1yXAY9RrRcSiKKRUK3D4lG3sq7Pj4YJ3oOEQU4VhUYkCT04OH39mPC7ItmJ7HnUMp9MZnmjE5y4Jf/+sQXF6f6DhEFMFYVKKcJElY8s4BdPkD+NFlI6BQcAIthZ5CocC86TlodHrwwudcrkxEg8eiEuXe3FWJDYfrce9lI2CJ14qOQzEk1aTHdydm4E+bS1Ha6BQdh4giFItKFDte347H/nkQs0ZacXFOkug4FIO+OzETyQYtlqw7gECAE2uJaOBYVKKU2+vHgtf2wGrQYd70HNFxKEZp1Urce+kIbC9rwV+/PiE6DhFFIBaVKLX0H9+gosWFB64qgE7NbfJJnHGZZlwzJhVPf3SEl4CIaMBYVKLQ2u0VeGt3FebPyMWwxHjRcYhw+5RsWBK0+PmbRfD6AqLjEFEEYVGJMnsqWrH0H9/g6tGpmFnIe/mQPOg1Kiy8Ig8Hqx1Y/tFh0XGIKIKwqESRhvZO/PRvu5FnNWDetOGi4xD1kW8z4oeXDMfqr07gX/tqRMchogjBohIlOrv8+NGru+D1BfDAVQVQq/hHS/JzzZhUzMhPxkNv78c31W2i4xBRBOCnWRSQJAm/fHs/Dte043+uGYmkBO6XQvKkUChw76UjkGHRY97qHShv7hAdiYhkjkUlCjz/WTH+sa8G983MQ57VIDoO0VnpNSo8dO0oaFVK3PXKDjS0d4qOREQyxqIS4d7eXYUVG47h1guHYRrv40MRwhSnwZLrR8Hp8eHWl75GZYtLdCQikikWlQj25fEm/PKd/bhylA03T84UHYdoQKxGPR67cQw8vgBuWbUVh2sdoiMRkQyxqESooko7fvTqLozPNGH+jBzebJAiks3UXVYMOjVufuErvLGjApLErfaJ6N9YVCLQ8fp2zPvLDmQlxeFnVxVCreQfI0UuS7wWj80Zg+n5KXh43QEseG0Pauxu0bGISCb4CRdhypo6cOf/2w5znAa/uHYU9Bpuj0+RT6dW4UeXjcADVxZga0kzZv1+E55efwQNDk60JYp1CimCx1kdDgfMZjPa2tpgMplExwm5E00duO3lr6FWKvHIDaNhiecyZIo+bq8f7++vwQcHauHzS7hytA3XjU3DJXnJyLTEiY5HREEwkM9vFpUIUdroxO1/3gbVyZKSyJJCUa7D48NXJU3YfKwRpY3d+60kJWiRlRSHTEsc4rVq6DVKSBLgD0jw+ALo7PLD4wvA4/PD55egVCigVilgitMgMV6DdHMcRqQkoCDViDxrAud2EQnCohJlvqluw9y/7EC8VoVfzWZJodjT3tmFI7XtKG9xobG9E62uLnh8fnT5JSgAKJWARqWEVqWERqWEWqWASqFAAN0lxuX1ocPjR0N7Jzo8fgCAOU6Di3ISceUoG74zOhU2k17o75EolrCoRJGtxU340au7kG6Jw0PXjoRRrxEdiShiSZKE9k4fTjR34Fh9Ow7VOnC0rh0BCbhkRDJuvXAYrh+fhnitWnRUoqjGohIlXttejqX/OIhxGSYsurqQE2eJQqC9swt7KlrxxfEmHKxxwKhX4/Yp2bjrkuHISooXHY8oKrGoRDiPz48nPjiMV78uxzVjUjF3Wg5USl5LJwq1ekcnPj1cj8+PNsDt9eOG8em474o8jM0wi45GFFVYVCLYiaYOLFy7B8fq23HXJcPxnTFpoiMRxZzOLj82H2vEhwdq0dDuwcxCK+6/Mh8X5ySJjkYUFVhUIlAgIOFv28vx1EdHYIrT4IErC5CbkiA6FlFM8wckfF3ajH8WVaOy1Y2LcxKxYFY+rii0csUQ0RCwqESYgzVtWPqPg9hd3oqrRtlwx9RsTuYjkpGAJGFPeSv+ua8GxxucGJlmxI8vG4E5EzOgVXPfTKKBYlGJEJUtLqzYcAzv7a1GhiUO91yaizHpkff7IIoVkiThUK0DH+yvxd5KO6xGHe6Yko3bp2QjzRw7y5sDJ/etkSBBAQW0aiXn0Z2DJHWfs4DUfc40KgXUqtgtuRFXVF544QX83//9H+rq6jBx4kQ8//zzmDJlyjm/LxKLiiRJ2FtpxytflOGjb2ph0mvwHxcOw6yRNv6HThRBqlpd+PhgHb4qbobXF8BlhSn4jwuG4TtjUiN6hV6bqwslTU6UN3egvNmF6lY3ats6Ue/ohN3dBYe7Cx5f4LTv06uVMOjVSE7QwWrUIcOiR6YlHtnJcchJTsCIFAPM8dG5vYLL60NpYwdOnDxnVa1u1LW5UefoRGtHF+wuLzr7OWdalRIGnQpJCTokG7TItMQhwxKH7KR4DE+OxwirASkGbVReZoyoovL3v/8dc+fOxUsvvYSpU6di5cqVeOutt3D06FHYbLazfm+kFBVJknC8wYkNh+rxzp4qlDZ2IM2kx/Xj0zCz0AqdOnL/UiOKdS6vD18VN+PL4kYcq3dCr1HiikIrrhqdissKrLIdaWnt8OJYfTuONThxvL4dx+qdON7Qjmant/cYS7wGKQYdkhK0sMRpYNSrYdCpoVOroFEroUD3ZTFfQEJnlx8dHj/aO7tgd3ehpcOLpnYP7O6u3tdLMWiRbzNiZKoBBalGFKYaUZhqiJjbgTg9PhQ3OHGsrr373J08f3Vt/74nlUGnhtXYfc4S4zUw6jUw6NTQa1TQnjxnANDlD8DjC8Dl9cPR2YU2VxdaXF40Oz1odnrR88Fs0qtPnisDCmz/PmdWoy6iC0xEFZWpU6fi4osvxh//+EcAQCAQQFZWFv77v/8bDz/88Fm/V65FxePzo6ShA/ur7NhV3oqvS5pRbXdDp1bi4pwkXFaQgnEZZig5gkIUVeraOrGjrBk7y1tR3OAEAOSmJODC4YmYlGXB6HQTClMNYdu40enxoarVhYpmFypaXChp7EBpoxPHG5xo6eguJCqlAhlmPTIscchM7L49Qbo5DmkmPeK0Q/9HlNvrR52jEzV2N2rsblSd/N8auxuBk58+yQYt8q0G5NsMyE1JQE5yArKS4jEsMQ4JuvDO1+vs8qPa7kZFS/d5K2vqQMnJc9ZTSBQAbCY9hp1yzjIseqSZ42AIQl6vL4B6Rydq2zpRbXejutXV/b92N7r83SfNpFcjz2ZAvtWAXGsCck+es6ykeJjj5D9yFTFFxev1Ij4+Hm+//Ta+973v9T4/b9482O12/OMf/+hzvMfjgcfj6f26ra0N2dnZqKysDElRkSQJXX4JvkAAnq4AvL4AXF0+uL1+tHf64Ojsgt3VhWanF43OTlS1ursfLW74JQlKBZCVFId8mxHjM80YlWaCNoavSRLFEsfJbf+P1DtwoqkDFS2u3g/mpAQNMi3xyLTokWLUI9mggSlOA5NOg3idCnp197++1UpF7yVhSQK6AgH4fBI8/gDcXX64vT44PT443D7Y3d2XGJravWho70R9eyecnf7ePFq1EulmPWxGHdLNccgw93y46qFWhv/vJV8ggLqTH8S1bf++vFTn8MB7ymUSg16FVKMeNqMeyQYtkhK0MOk1MMerkaBVI+7kPZ90KhU0akClVKLn34CS1P0+Pj/g8XffC6rT64fT64PD3YU2d8/f4R40Oj2ob+uE3e3rfW+1SoFUow42ox7pJ8tcuqX7f3Wq8I+EByQJDe2dqLF3F7/atk7UOdyod3TC5f33OYvXKpFmikOqqfucJRu0sMRrYInTdo+K6TVI0CoRp1FDr1VBp1JCp1ZCoz55C4qTP3ehHLFxOBzIysqC3W6H2Xz2fYqELi1pamqC3+9Hampqn+dTU1Nx5MiR045fvnw5Hn/88dOez8rKClnGoSoH8KXoEEQkK5UA9gl43xIB7xkMhwW+d5nA9x6Ko6IDnKf29nZ5F5WBWrJkCRYvXtz7dSAQQEtLC5KTkyP6Wt2pelpmqEaJIg3PR188H6fjOemL56Mvno++5HI+JElCe3s7MjIyznms0KKSkpIClUqF+vr6Ps/X19cjLe30HVl1Oh10Ol2f5ywWSygjCmMymfgf1Sl4Pvri+Tgdz0lfPB998Xz0JYfzca6RlB5CJ0xotVpceOGF2LhxY+9zgUAAGzduxLRp0wQmIyIiIjkQfuln8eLFmDdvHi666CJMmTIFK1euREdHB+bPny86GhEREQkmvKjcdtttaGxsxNKlS1FXV4dJkyZh/fr1p02wjRU6nQ6PPfbYaZe4YhXPR188H6fjOemL56Mvno++IvF8CN9HhYiIiOhMuKkHERERyRaLChEREckWiwoRERHJFosKERERyRaLigAvvPACcnJyoNfrMXXqVOzYseOsx69cuRIjR45EXFwcsrKy8POf/xydnZ1n/Z5IMpDz0dXVhV//+tfIy8uDXq/HxIkTsX79+jCmDa0tW7Zgzpw5yMjIgEKhwHvvvXfO79m0aRMuuOAC6HQ65OfnY82aNSHPGS4DPR+1tbW44447UFhYCKVSiUWLFoUlZ7gM9HysW7cO3/nOd2C1WmEymTBt2jR8/PHH4QkbBgM9H19++SVmzJiB5ORkxMXFYdSoUfjDH/4QnrBhMJi/P3p89dVXUKvVmDRpUsjyDRaLSpj9/e9/x+LFi/HYY49hz549mDhxIq699lo0NDT0e/zatWvx8MMP47HHHsPhw4fxyiuv4O9//zt+9atfhTl5aAz0fDzyyCP405/+hOeffx6HDh3Cfffdh5tvvhl79+4Nc/LQ6OjowMSJE/HCCy+c1/FlZWW44YYbMGvWLBQVFWHRokW49957o+bDaKDnw+PxwGq14pFHHsHEiRNDnC78Bno+tmzZgu985zv48MMPsXv3bsyaNQtz5syJ2f9eEhIScP/992PLli04fPgwHnnkETzyyCN4+eWXQ5w0PAZ6PnrY7XbMnTsXV111VYiSDZFEYTVlyhRp4cKFvV/7/X4pIyNDWr58eb/HL1y4ULryyiv7PLd48WJpxowZIc0ZLgM9H+np6dIf//jHPs/dcsst0p133hnSnCIAkN59992zHvPQQw9JY8eO7fPcbbfdJl177bUhTCbG+ZyPU82cOVP62c9+FrI8og30fPQYM2aM9Pjjjwc/kGCDPR8333yz9MMf/jD4gQQbyPm47bbbpEceeUR67LHHpIkTJ4Y012BwRCWMvF4vdu/ejauvvrr3OaVSiauvvhpff/11v98zffp07N69u/dySGlpKT788EPMnj07LJlDaTDnw+PxQK/X93kuLi4OX34Zm/eo/vrrr/ucPwC49tprz3j+KLYFAgG0t7cjKSlJdBRZ2Lt3L7Zu3YqZM2eKjiLM6tWrUVpaiscee0x0lDMSvjNtLGlqaoLf7z9t193U1FQcOXKk3++544470NTUhEsvvRSSJMHn8+G+++6Liks/gzkf1157LVasWIHLL78ceXl52LhxI9atWwe/3x+OyLJTV1fX7/lzOBxwu92Ii4sTlIzk6Pe//z2cTid+8IMfiI4i1LBhw9DY2Aifz4dly5bh3nvvFR1JiOPHj+Phhx/GF198AbVavnWAIyoyt2nTJjz55JN48cUXsWfPHqxbtw4ffPABfvOb34iOJsSzzz6LgoICjBo1ClqtFvfffz/mz58PpZI/ykRns3btWjz++ON48803YbPZRMcR6osvvsCuXbvw0ksvYeXKlXj99ddFRwo7v9+PO+64A48//jgKCwtFxzkr+VaoKJSSkgKVSoX6+vo+z9fX1yMtLa3f73n00Udx11139Tb+8ePHo6OjAz/+8Y/xv//7vxH9AT2Y82G1WvHee++hs7MTzc3NyMjIwMMPP4wRI0aEI7LspKWl9Xv+TCYTR1Oo1xtvvIF7770Xb7311mmXCmNRbm4ugO6/T+vr67Fs2TLcfvvtglOFV3t7O3bt2oW9e/fi/vvvB9B9aVCSJKjVanzyySe48sorBafsFrmfchFIq9XiwgsvxMaNG3ufCwQC2LhxI6ZNm9bv97hcrtPKiEqlAgBIEX6bpsGcjx56vR6ZmZnw+Xx45513cNNNN4U6rixNmzatz/kDgA0bNpzz/FHseP311zF//ny8/vrruOGGG0THkZ1AIACPxyM6RtiZTCYcOHAARUVFvY/77rsPI0eORFFREaZOnSo6Yi+OqITZ4sWLMW/ePFx00UWYMmUKVq5ciY6ODsyfPx8AMHfuXGRmZmL58uUAgDlz5mDFihWYPHkypk6diuLiYjz66KOYM2dOb2GJZAM9H9u3b0d1dTUmTZqE6upqLFu2DIFAAA899JDI30bQOJ1OFBcX935dVlaGoqIiJCUlITs7G0uWLEF1dTVeffVVAMB9992HP/7xj3jooYdwzz334LPPPsObb76JDz74QNRvIagGej4AoKioqPd7GxsbUVRUBK1WizFjxoQ7ftAN9HysXbsW8+bNw7PPPoupU6eirq4OQPcEdLPZLOT3EEwDPR8vvPACsrOzMWrUKADdy7d///vf44EHHhCSP9gGcj6USiXGjRvX5/ttNhv0ev1pzwsneNVRTHr++eel7OxsSavVSlOmTJG2bdvW+2szZ86U5s2b1/t1V1eXtGzZMikvL0/S6/VSVlaWtGDBAqm1tTX8wUNkIOdj06ZN0ujRoyWdTiclJydLd911l1RdXS0gdWh8/vnnEoDTHj3nYN68edLMmTNP+55JkyZJWq1WGjFihLR69eqw5w6VwZyP/o4fPnx42LOHwkDPx8yZM896fKQb6Pl47rnnpLFjx0rx8fGSyWSSJk+eLL344ouS3+8X8xsIssH893IquS5PVkhShF8/ICIioqjFOSpEREQkWywqREREJFssKkRERCRbLCpEREQkWywqREREJFssKkRERCRbLCpEREQkWywqREREJFssKkRERCRbLCpEREQkWywqREREJFssKkRERCRb/z+QFzc6JQkYEgAAAABJRU5ErkJggg==",
      "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": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}