Synthetic Control Simulation V1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "70e24549",
   "metadata": {},
   "source": [
    "# Synthetic Control Simulation\n",
    "\n",
    "Sources\n",
    "1. Causal Inference Mixtape\n",
    "2. Causal Inferece for the Brave and True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "069dd452",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "np.random.seed(777)\n",
    "\n",
    "from toolz import partial\n",
    "from scipy.optimize import fmin_slsqp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "431895d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/alexander.clark\n"
     ]
    }
   ],
   "source": [
    "%cd ~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "eff3da1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from synth_control import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1bbf593",
   "metadata": {},
   "source": [
    "## Generate Fake Data\n",
    "\n",
    "For each unit $i$ and time period $t$, we observe features $A_{it}, B_{it},$ and $Y_{it}$. $Y_{it}$ is the outcome of interest and formed as a linear combination of $A_{it}$ and $B_{it}$. $A_{it}$ is drawn randomly at $t=0$ and $A_{it} = A_{it-1} + \\epsilon_{it}$ for $t\\geq 1$, where $\\epsilon_{it}$ is normal noise. The same is true of $B_{it}$, with noise iid across units, time, and features.\n",
    "\n",
    "The treated unit is a convex combo (`true_weights`) of the donor units, plus some random noise. \n",
    "\n",
    "In treatment periods, $Y_{it}$ is shifted up by 0.5 (`treatment_effect`) for the treated unit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "55764c19",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = 'ABY'\n",
    "donor_units = range(3)\n",
    "treated_unit = 'T'\n",
    "\n",
    "n_periods = 10\n",
    "periods = range(n_periods)\n",
    "n_treament_periods = 3\n",
    "pretreatment_periods = range(0, n_periods - n_treament_periods)\n",
    "treatment_periods = range(n_periods - n_treament_periods, n_periods) # or post-treatment\n",
    "\n",
    "treatment_effect = 0.5 # additive effect on Y "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a5c3e955",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.DataFrame()\n",
    "\n",
    "\n",
    "# Construct Dataset for the donor units\n",
    "for i in donor_units:\n",
    "    for feature in features[0:-1]: # exclude the outcome feature\n",
    "        for period in periods:\n",
    "            \n",
    "            # initialize at randomly\n",
    "            if period == 0:\n",
    "                data.loc[i,feature + str(period)] = np.random.normal() \n",
    "            \n",
    "            # evolve according to random shock\n",
    "            else:\n",
    "                prev = data.loc[i,feature + str(period-1)]\n",
    "                data.loc[i,feature + str(period)] = prev + np.random.normal(0,1)\n",
    "      \n",
    "    # Generate outcome variable\n",
    "    for period in periods:\n",
    "        # outcome is a linear combo of a and b\n",
    "        data.loc[i,'Y' + str(period)] = .1*data.loc[i,'A'+str(period)] \\\n",
    "                            + .9*data.loc[i,'B'+str(period)] + np.random.normal(0, .01) # adjust var\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a6fe4956",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the Treated Object\n",
    "\n",
    "# Untreated potential outcome as a combo of donors\n",
    "true_weights = np.random.dirichlet(np.ones(len(donor_units)))\n",
    "\n",
    "m = np.matrix(data)\n",
    "convex_combo = np.array(np.dot(true_weights, m)).flatten()\n",
    "\n",
    "# Think of this as Y_0 | D = 1, the potential outcome of no treatment for the treated unit\n",
    "data.loc[treated_unit] = convex_combo + np.random.normal(0, .01, size = len(convex_combo)) # adjust var\n",
    "\n",
    "# Treated outcome for treated unit\n",
    "data.loc[treated_unit + 'treat'] = data.loc[treated_unit]\n",
    "for treatment_period in treatment_periods:\n",
    "    data.loc[treated_unit + 'treat','Y'+str(treatment_period)] += treatment_effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9aa1ab87",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get dataframe with just pretreatment periods \n",
    "\n",
    "pretreatment_cols = []\n",
    "\n",
    "for i in list(data):\n",
    "    i2 = i\n",
    "    for x in features:\n",
    "        i2 = i2.lstrip(x)\n",
    "    if int(i2) in pretreatment_periods:\n",
    "        pretreatment_cols.append(i)\n",
    "        \n",
    "pretreatment_data = data[pretreatment_cols]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c820f39a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/HElEQVR4nO3dd3wU1frH8c+z6YGQQAgt9BpagBB6R1EBBUURexdFBOu1/a5XvNbr1WtXxIqKAiI2QLoIKi2hQ+g11JBAIJC+5/fHLhgxgSS7m9lsnvfrtS+yu7Nzvjtsnp2cM3NGjDEopZTyfTarAyillCobWvCVUqqC0IKvlFIVhBZ8pZSqILTgK6VUBaEFXymlKggt+MqriEgLEVktIidFZKyIhIjITyKSLiLfWJ3Pl4jIOBH50oPrzxCRxp5avyo5LfgVlPOX8czNLiKZBe7f6MF2bxOR/HPazxCROs5FHgMWGWPCjDFvAdcANYFIY8xwF9r1aHEr0E6vAu/plIiYc95nfTe2tVtELnbX+s5Z920i8psrbRpjKhtjdjpf95mIPO/unKpk/K0OoKxhjKl85mcR2Q3cZYyZf+5yIuJvjMlzc/NLjTE9i3iuATD5nPtbPZDBI4wxS4DKACLSENgFRBSW30PbVqki6R6++gsR6SsiySLyuIgcAj4tbG/Puefa1PlzkIi8KiJ7ReSwiIwXkZBStL0Q6Ae849wb/hr4FzDCef9O53J3iEiSiBwTkTki0qDAOlqLyDwRSXNmeUpELgOeKrCetUW031JEFonIcRHZKCJDCjz3mYi8KyIznd1Ny0WkSQnf3zgRmSYiX4rICeA2EQkXkY9F5KCI7BeR50XEz7l8ExFZKCKpInJURCaJSITzuS+A+sBPzvf0mPPxriLyh/M9rBWRvgXabyQivzrzzwOqlyR/Ie/nvNvkzGdEREYCNwKPObP+5Eq7qvS04KvC1AKq4di7HlmM5f8DNAfaA02BaByFukSMMf2BJcD9zu6A64EXgSnO+x+LyJU4ivcwIMq5/NcAIhIGzAdmA3WcWRYYY2afs55257YtIgHAT8BcoAYwBpgkIi0KLHY98CxQFdgOvFDS9wgMBaYBEcAkYCKQ58zaAbgEuOtMLOAl53tpCdQDxjm31c3AXuAK53t6RUSigZnA8zj+/x4FvhWRKOf6vgIScRT654BbS5H/XBfcJsaYCc73+ooz6xVuaFeVghZ8VRg78IwxJtsYk3m+BUVEgLuBh4wxacaYkziK63XneVlX5x7omduOEmS7B3jJGJPk7A55EWjv3Mu/HDhkjHnNGJNljDlpjFlezPV2xdEV87IxJscYsxCYgaOgnTHdGLPC2e4kHF9wJbXUGPO9McYOVAEGAg8aY04ZY44Ar+PcdsaY7caYec7/hxTgf0Cf86z7JmCWMWaWMcZujJkHJACDnGMHnYCnnetbjOMLzlXu2CaqjGgfvipMijEmq5jLRgGhQKKj9gOOPVO/87xm2Xn68C+kAfCmiLxW4DHB8VdFPaAkXx4F1QH2OQvxGXuc6z3jUIGfT+Psqy+hfQV+bgAEAAcLbDvbmWVEpAbwFtALCHM+d+w8624ADBeRgnvQAcAvON7fMWPMqQLP7cGxzQqT53ztuQKA3AL33bFNVBnRgq8Kc+4UqqdwFHUARKRWgeeOAplAa2PM/jLItg94wRgz6dwnnHv51//9JcDf39O5DgD1RMRWoOjXB7aWOumFc+wDsoHqRQzevuRcPtYYk+rsznqniHWdWd8Xxpi7z12Rc9tUFZFKBYp+/ULWccZeoL6IiHFOqSsioTi6u/ac7w0WQafl9QLapaOKYy3QWkTai0gwzn5kAGdx/BB43blHiohEi8ilHsoyHnhSRFo72woXkTOHa84AaonIg86B5DAR6eJ87jDQUESK+swvx/HF9piIBDgHO6/gr0cMuZUx5iCOMYPXRKSKiNicA7Vnum3CgAzguLN//h/nrOIwUPA49y+BK0TkUhHxE5FgcQzC1zXG7MHRvfOsiASKSE/n+yvKciALeMK5nkrAy851lKbgn5tVWUALvrogY8xW4N84BkS3Aecen/04jgG7Zc6jT+YDLShaN/n7cfidipnlOxyDxJOdbW3A0Q+Oc/xgAI5CdsiZtZ/zpWdO2koVkVWFrDcHGOJc11HgPeAWY8zm4uRywS1AILAJR3fNNKC287lngTggHcdg7PRzXvsS8E/nOMijxph9OAaFnwJScOzx/4M/f89vALoAacAzwOdFhTLGZAODgb5AMrATR7fQtWf2+EvoY6CVM+v3pXi9cgPRC6AopVTFoHv4SilVQWjBV0qpCkILvlJKVRBa8JVSqoJw6Th8EakGTAEaArtxjOD/7cQQ5/wfHwFtcByPe4cxZumF1l+9enXTsGFDVyIqpVSFkpiYeNQYE1XYc66eePUEjrlKXhaRJ5z3Hy9kuTeB2caYa0QkkAIn8ZxPw4YNSUhIcDGiUkpVHCJS5HkSrnbpDMUx+RPOf68spPEqQG8cx+HinKfkuIvtKqWUKiFXC35N59mCZ84arFHIMo1xnATyqTiuZPSR86y9QonISBFJEJGElJQUF+MppZQ644IFX0Tmi8iGQm5Di9mGP46zBd83xnTAcfr6E0UtbIyZYIyJN8bER0UV2g2llFKqFC7Yh2+MKfJyZuK4wERtY8xBEakNHClksWQgucA0tdM4T8FXSvme3NxckpOTycoq7iSs6kKCg4OpW7cuAQGFTWpaOFcHbX/EcRGFl53//nDuAsaYQyKyT0RaGGO2ABfhmDdEKVVBJCcnExYWRsOGDSkwFbQqJWMMqampJCcn06hRo2K/ztU+/JeBASKyDcekVS8DiEgdEZlVYLkzVw9ah+MCCS+62K5SqhzJysoiMjJSi72biAiRkZEl/ovJpT18Y0wqjj32cx8/AAwqcH8NEO9KW0qp8k2LvXuVZnv63pm2uZnwx9vYdy62OolSSnkV3yv4Nn9O/vImm7551uokSikvsnv3btq0afOXx8aNG8err75a5GvWrFnDrFmziny+pEq7vr59+7rlJFTfK/h+AWyoczVtMhPYs2WN1WmUUuXY+Qp0Xl5hV6Us/frKgu8VfKDF4LFkmwAOzXvT6ihKqXKgb9++PP7443Tu3JnmzZuzZMkScnJy+Ne//sWUKVNo3749U6ZMYdy4cYwcOZJLLrmEW265hZSUFK6++mo6depEp06d+P333wFYsWIF3bt3p0OHDnTv3p0tW7YUur5Tp05xxx130KlTJzp06MAPPzgOdMzMzOS6664jNjaWESNGkJmZ6Zb36ZMXMa9WI5qVVS+mbcpM0o8dJbxqdasjKaUKGPHB3+dOvDy2Njd3a0hmTj63fbrib89f07Euw+PrkXYqh1FfJv7luSn3dHM5U15eHitWrGDWrFk8++yzzJ8/n3//+98kJCTwzjuOa8ePGzeOxMREfvvtN0JCQrjhhht46KGH6NmzJ3v37uXSSy8lKSmJmJgYFi9ejL+/P/Pnz+epp57i22+//dv6nnrqKfr3788nn3zC8ePH6dy5MxdffDEffPABoaGhrFu3jnXr1hEXF+fy+wMfLfgAVfuNIfS7n1k26z263vgvq+MopSxW1FEtZx4fNmwYAB07dmT37t1FrmfIkCGEhIQAMH/+fDZt+vO0ohMnTnDy5EnS09O59dZb2bZtGyJCbm5uoeuaO3cuP/7449lxhKysLPbu3cvixYsZO3YsALGxscTGxpbszRbBZwt+03Y9OPRLHPGHvwH7/4HNz+pISimn8+2RhwT6nff5apUCS7VHHxkZybFjf529PS0t7eyJS0FBQQD4+fmdt3++UqU/pwKz2+0sXbr07BfAGWPGjKFfv35899137N69m759+xa6LmMM3377LS1atPjbc544jNUn+/DPqDXgAfxP7IWtc6yOopSyWOXKlalduzYLFiwAHMV+9uzZ9OzZs8jXhIWFcfLkySKfv+SSS852z4BjUBYgPT2d6OhoAD777LMi13fppZfy9ttvY4wBYPXq1QD07t2bSZMmAbBhwwbWrVtXgndaNJ8u+MRcTnalOuz5+TWrkyilvMDnn3/O888/T/v27enfvz/PPPMMTZo0KXL5fv36sWnTprODrOd66623SEhIIDY2llatWjF+/HgAHnvsMZ588kl69OhBfn5+ket7+umnyc3NJTY2ljZt2vD0008DMGrUKDIyMoiNjeWVV16hc+fObnn/cuabxRvFx8cbV489Xfr503Tb+Ra7rp1Po1ad3JRMKVUSSUlJtGzZ0uoYPqew7SoiicaYQmc28O09fKDloNFkmkCOzH/L6ihKKWUpny/4EdVrsT7yUmJTZ3P86GGr4yillGV8vuADRF00hhDJIWnWOxdeWCmlfFSFKPiNWndha2gcsQemQn7JT4dWSilfUCEKPkDzIY9SKesQbJlpdRSllLJEhSn4NL8ME9GA9EVvW51EKaUsUXEKvs2PhJrXEH5kJdvX/WF1GqVUGUpNTaV9+/a0b9+eWrVqER0dffZ+Tk5Oqda5aNEi/vij5LWkYcOGHD16tFRtuspnp1YoTPPL7uPU5ndIW/g2xHa3Oo5SqoxERkaePQt23LhxVK5cmUcfffTs83l5efj7l6wcLlq0iMqVK9O9e/mpJRWq4IdXrc6yqEF0SJlB2pH9VKsRbXUkpZRFbrvtNqpVq8bq1auJi4vjvvvuY/To0aSkpBAaGsqHH35ITEwMP/30E88//zw5OTlERkYyadIkMjMzGT9+PH5+fnz55Ze8/fbbxMTEcO+997J3714A3njjDXr06EFqairXX389KSkpdO7cGStPdnWp4ItINWAK0BDYDVxrjDlWyHIPAXcBBlgP3G6MKdnVd92k9sVjCZr8HatmvU232162IoJSFdvPT8Ch9e5dZ622MLDkv89bt25l/vz5+Pn5cdFFFzF+/HiaNWvG8uXLue+++1i4cCE9e/Zk2bJliAgfffQRr7zyCq+99hr33nvvX/5SKGqq5GeffZaePXvyr3/9i5kzZzJhwgT3vvcScHUP/wlggTHmZRF5wnn/8YILiEg0MBZoZYzJFJGpwHXAZy62XSoNYuJYF9yRFnunQv5z4BdgRQyllBcYPnw4fn5+ZGRk8McffzB8+PCzz2VnZwOQnJzMiBEjOHjwIDk5OWdn1zxXUVMlL168mOnTpwMwePBgqlat6sF3dH6uFvyhQF/nzxOBRZxT8Au0EyIiuUAocMDFdl3SaPAjhH17A2z6AdpeY2UUpSqeUuyJe8qZqY7tdjsRERFn+/kLGjNmDA8//DBDhgxh0aJFjBs3rtB1FTVVMnhmquPScPUonZrGmIMAzn9rnLuAMWY/8CqwFzgIpBtj5ha1QhEZKSIJIpKQkpLiYrzChbUeCNUaY1823iPrV0qVL1WqVKFRo0Z88803gGOe+rVr1wJ/nep44sSJZ19z7lTHRU2VXHCq459//vlvc/KXpQsWfBGZLyIbCrkNLU4DIlIVx18CjYA6QCURuamo5Y0xE4wx8caY+KioqOK+j5Kx2djT9CZs+1eyddWvnmlDKVWuTJo0iY8//ph27drRunXrs9eXHTduHMOHD6dXr15Ur/7n5VKvuOIKvvvuO9q3b8+SJUuKnCr5mWeeYfHixcTFxTF37lzq169vyfsDF6dHFpEtQF9jzEERqQ0sMsa0OGeZ4cBlxpg7nfdvAboaY+670PrdMT1yUU6mpyH/a8Xm8J7EPzzNI20opRx0emTPKOvpkX8EbnX+fCvwQyHL7AW6ikioODqyLgKSXGzXZWHh1dhQ43Ji0xdy9OBeq+MopZTHuVrwXwYGiMg2YIDzPiJSR0RmARhjlgPTgFU4Dsm0AdYdl1RA9KUP4o+dbbN0rnyllO9z6SgdY0wqjj32cx8/AAwqcP8Z4BlX2vKEek3bsDa0M832fUNO1nMEBv99dF0p5R7GGK85WsUXlKY7vuLMpVOE8H5jqM5xArcU1hullHKH4OBgUlNTLT3L1JcYY0hNTSU4OLhEr6tQUysUpmGny2FFC1j2PsSOAN0DUcrt6tatS3JyMp461LoiCg4Opm7duiV6TYUv+IiQ0f4OKs9/nC2JC2gRf7HViZTyOQEBAUWeoarKToXv0gGQdtdzglAyftVLICqlfJcWfKBSWDibag6l3YlfOZy8w+o4SinlEVrwnepf9iCCYefPeoimUso3acF3qtMohnWVuhOz/1uyTmdYHUcppdxOC34BwT1HU5WT2NfrVAtKKd+jBb+Alt0GQY3WhK76CPR4YaWUj9GCX5AIpss9cHg9u1fNszqNUkq5lRb8c2S3HMZxKpO2UAdvlVK+RQv+OYJDw0iqPYx2Gb9xcM8Wq+MopZTbaMEvRKNBD2AQds/WvXyllO/Qgl+IWvWasjasF60OfsfpjHSr4yillFtowS9CaK/RhHOKtGWTrI6ilFJuoQW/CDGdBmCvFUvdLZ/rIZpKKZ+gBb8IYrNh6zoKUpJI3zTf6jhKKeUyLfjnYVpfxXEJZ+fM16yOopRSLnOp4IvIcBHZKCJ2ESn0KunO5S4TkS0isl1EnnClzbIkASEkRV9Du1PL2L9zo9VxlFLKJa7u4W8AhgGLi1pARPyAd4GBQCvgehFp5WK7ZabpoLHkY2PfHD1EUylVvrlU8I0xScaYC52d1BnYbozZaYzJASYDQ11ptyxF1WnI2ip9aX3oB06dOGZ1HKWUKrWy6MOPBvYVuJ/sfKxQIjJSRBJEJMFbrn8Z1ud+wiSTzXMmWB1FKaVK7YIFX0Tmi8iGQm7F3Usv7KrgRR7naIyZYIyJN8bER0VFFbMJz2oR35+Tke2IO/QN2O1Wx1FKqVK54EXMjTGuXtU7GahX4H5d4ICL6yxzYX3uh+l3w46F0EwvdK6UKn/KoktnJdBMRBqJSCBwHfBjGbTrXq2u5HRgdZK+f8XqJEopVSquHpZ5lYgkA92AmSIyx/l4HRGZBWCMyQPuB+YAScBUY0z5O8bRP5C1ta6m5anl7N26xuo0SilVYq4epfOdMaauMSbIGFPTGHOp8/EDxphBBZabZYxpboxpYox5wdXQVmk2aAw5xp+Dc/UQTaVU+aNn2pZA9Vr1WBtxEW1SZnLieKrVcZRSqkS04JdQRL8xVJIsNs16z+ooSilVIlrwS6hZ+17sD2tHh4NTwZ5vdRyllCo2LfilEH3pgwSd3Avb5lodRSmlik0Lfmm0vIK8yrXZP+d1q5MopVSxacEvDb8AVkYNIzptObuTEq1Oo5RSxaIFv5RiBo0hywRweP6bVkdRSqli0YJfSlWjarOu2iXEHv2Z9DTvmORNKaXORwu+C6pf9AAhkkPSzLetjqKUUhekBd8Fjdt0ISmoHa2Sp0B+ntVxlFLqvLTguyhm6KNUyT4EW3+2OopSSp2XFnwXScxgTHg9Ti9+x+ooSil1XlrwXWXzY12dawk9uIzt65danUYppYqkBd8NGg4YxWkTRNoCnUVTKeW9tOC7QXi1KNZXH0i7Y/M4lnLQ6jhKKVUoLfhuUuvisQRJLptn6SGaSinvpAXfTRq07MiGoA403T0Zk5djdRyvk3HiGLs2rbQ6hvJGRzbD+z1g12Krk/g8LfhuVOuSh4gyqcjmGVZH8SrpqYc5/GZ/6k8ZwKZls62Oo7zNgn/D4Q3wzW1wfK/VaXyaFnw3qt7hCqjaCLN8vNVRvEfmcdInXE69vL0clWpUmz2a9LQjVqdS3iI5EbbMhA43OU5enHwj5GZancpnuXoR8+EislFE7CISX8Qy9UTkFxFJci77gCttejWbjUMxtyD7lrN9tf55StYJ+PJq6uXuYlu/8aQP+ZRIc4xDX44EY6xOp7zBwucgNBIuexmu/hAOrYefHtDPh4e4uoe/ARgGnK+65QGPGGNaAl2B0SLSysV2vVZol1s5ZYI5tqhin4h16uRxDr43GHNwDXLtRFr3HU7zuD4cjH+UFmm/QOJnVkdUVtv9G+z8BXo+BEFh0PxS6PcUrJsC+leyR7hU8I0xScaYLRdY5qAxZpXz55NAEhDtSrverEpEJBuiBtPu+AKOHtpndRxLnM5IZ8/blxOVvoGdfd6CmMFnn6s/+Alo3Bcz+0kObFttYUplKWNg4fNQuRZ0uuvPx3s9CjGXw5z/00FcDyjTPnwRaQh0AJafZ5mRIpIgIgkpKeVz2uE6lz5AoOSxrQIeopl1+iQ73x5Ci+wNrO38Ck363PjXBWw2zJXjOZEfQNbXt5GVecqaoMpa2xfA3qXQ+1EICPnzcZsNrnwfIps4B3Er5k6Tp1yw4IvIfBHZUMhtaEkaEpHKwLfAg8aYE0UtZ4yZYIyJN8bER0VFlaQJr1GvWTvWBnei2d4p5GRnWR2nzGRlnmLrW0NplbWW1R1fpOPguwtdTqrUZk+vV2ls383aT8aWcUplOWMcffcR9SHu1r8/H1wFrvsK8nNhig7iutMFC74x5mJjTJtCbj8UtxERCcBR7CcZY6a7Eri8qNR7NNU5jv+cxxwfXF+Xl03uVzcRm5VIYrtniR9y33kXj+0/gmU1rqVLyjTWLJhcRiGVV9g8Aw6ugT5PgH9g4ctUbwbDPoSDa2HGQzqI6yYe79IREQE+BpKMMf/zdHveomm3K6H7WGyrJpL6/mWkHk62OpLHmLxs+OY2wvYtJGPAq3QaVrwDsTrc8SY7/BrRYMk/OHpgt2dDKu9gz4eFL0BkM4gdcf5lW1wGfZ+CtV/D8g/KJp+Pc/WwzKtEJBnoBswUkTnOx+uIyCznYj2Am4H+IrLGeRvkUuryQAQueY5Tg9+nUspact/vw7Y1S6xO5Xa5Odmsef1q2DILBr1K5R6Fd+MUJig4FP/hn1JJsqky+36w2z2YVHmFDd9CShL0exL8/C+8fO9/QIvBMOcpx1E9yiWuHqXznTGmrjEmyBhT0xhzqfPxA8aYQc6ffzPGiDEm1hjT3nmbdf41+45KnW4gedj3ANT77ioSfvSdw83ycnNY9/YIOpxawrLmj0Ln4hf7MxrEdCDw8v8SuHcJ5ne9ILxPy8+FRS9BzbbQ6qrivcZmg6vGOwZxp96qg7gu0jNty0DTdj0JGPUrO4NiiF/1OMvHjyr3l0TMz8tjzds30PHkLyxr8gBdb3i69CuLu4WsZleQv+A5tq5a5LaMysus+QrSdkL//3MU8uI6M4iblw1TbvL9Qdz8XDi22yOr1oJfRiJr1qXZowtYXn0YXQ59BZOugdNpVscqFWPPZ9U7NxF/Yh5LG95H15v/7doKRcge9DpHqUboT/dwMr18bhd1HnnZ8OsrEN0Rml9W8tdXbwbDJjgGe315EDcr3VEbPhkI2RluX70W/DIUEBhEl/s/xVzxFuz5naz3+pS/GSTtdmTGQ3Q6/jNL691Nt9tecstqw6tGcXzgu9S2H2bLxyPdsk7lRRI/gxPJ0P9px/hWacQMgr5POgZxV0xwazyvkJ4Mn1zmGKvo/08Iquz2JrTgW0A63or9lhmcOpVBzSmDWT1notWRisXY7Zz87kFYNRF6PUK3O/7r1vXHdLmUFQ3uJv7EPBJ+fN+t61YWyjkFi1+Fhr2gcV/X1tX7MWgxCGY/6VuDuAfXwocXOYr+jdOgw40Xfk0paMG3iK1BF+x3/8K+gIZ0WDqWZR89jD0/3+pYRTJ2O8vH30PY+omkx41ybU/tPDrd/AKbAtrQatWz2I/udPv6lQVWTIBTRxx7ra5+Zs4M4lZr7BjETfeBw523zYNPB4HNH+6YDU36eawpLfgWiqrTkAaP/MKKiEF0Tf6Yda8NJivjmNWx/sbY7SyfcD9dj0xlWY1rqXL5ix4p9gD+AYHUuHUiIcGB2KbfCXoxmfItKx1+ewOaDoD6Xd2zzuDwcwZxy/HZ7AmfwlcjoFojuGs+1Gzt0ea04FssOKQSncZOYnnME7Q9vYKgzy6Bo9utjnWWsdtZ9tGDdD00yTHgfO8HSEmOsCiF6nWbIkPehgOr2DXtKY+2pTxs6XuQddxxZI47RTWHYR/AgdUw8+HyN4hrt8P8cTDjQWjSH27/GarU9nizWvC9gNhsdLnuSWy3fo+cTsU+oR9rF061OhYA26b+k24HJrKi2hV0GvWRx4v9Wa2Gsr7WVTRI+oj1i4s9i4fyJqfTYOm70PIKqNPB/euPGQx9Hoc1k2DlR+5fv6fkZcP0u+C316Hj7XD9ZMf00GVAC74XkUa9YeQiDkgN2v46kqUTn8JYefbp4v/SfPO77Kx7JfGjJ2Lz8yvT5pve/Db7/OpSe+FY0o7sL9O2lRv89jrkZEA/N+/dF9TnCWg+EGY/Abt/91w77nI6Db64ynHG8cXj4PLXi3fGsZtowfc2EfWJHLuI1VX60W3Xu6z+35Wczkgv8xhrJ49zzFceO4LGd3xS5sUeIKRSGHlXfUQVk8G+T2+39stPlczJQ7DiQ4i9Fmq09Fw7Npuja6dqQ/jmVkj34h2DtF3w8SWQvBKu/thx4RcPjYUVRQu+FwqpFEbcQ9+yrPFY2p9czKHXe5Oy77zXmXGrZZP+TbvNr7Op2gAY+h7Yyr7Yn9GkbVdWxTxCu8zlLJ/ysmU5VAkteQ3yc6DvE55v68wgbm6W9w7iJifCxwPgVArc/D20vcaSGFrwvZTYbHS95Tk29P2ImvajVP/qMtj5q8fbXT7lZbpue41VlXrTfNRXZfrnZlG6jHiCbRE96LLtdcc1T5V3O77XcfRJh5sch0+WhagWjsM1D6yCmY941yDu5pnw2WAICIU750HDHpZF0YLv5WL7XUPo6F+RyjUwX1zF6qkveqxrY/k3r9El6SVWh3an7QPT8A8oYq7yMiY2G83u/hwJrYaZdge5We4/5Vy50a+vOLoq+jxWtu22vNxxYtaaL71nEHfZeJh8I9Rs5TjsMqq5pXG04JcDUr0p3DWf3ZG96LDpPyS8eb3bLw2Yl/A5XTb+m7UhnWk19lsCAoPcun6XVaqO/coPMEe3sWrCKKvTqKKk7nBMkhZ/J4TXLfv2+z7pmKtn9hOw54+yb/8Me77jbODZjzvODL51BlSuYV0eJy345UVQGA1GTWdp/ZF0Sp/N3tf6cjh5h1tWbdZ8jf+MseQ27EuLsd8TFBzqlvW6m61pX5bXvpkuaT+y6udPrY6jCvPLi+AfBL0etqZ9mw2u+gAiGsDUW6wZxM057Wh72XvQ5V4Y8QUEesfvlBb8csTm50e3O/7L6u7vEp27F7+P+rM9YZ5L60yc8SHm+/uwN+xFwI2TCQ6p5Ka0nhF/+6ts9W9O0+X/x6G926yOowo6vNFxuGGXe6zdmw2JcA7iZsLUm8t2EDcjBSZe4ei3v+xlGPgfSw96OJcW/HKowyU3cfS6meT6hdJk1vWOAbJSWPXzp7Rb+RibA1uTfc0kCAhxc1L3CwgMotINn2EzdtK+uI28XJ16wWv88qLjBKLuXnBh+hoxjkHc/Ykwq4wGcY9uh48vhsMbHHv1Xb2v61ELfjnVoGU8tR75w3Gy1owH2frRneRkF39PZvXcL2m77BG2BcbQYMwMQipX8WBa94pu3JrNHZ+hVe4G0ue6Z3pm5aL9iY6Lk3e7H0KrWZ3GoeUVjkskrv4SEj72bFt7ljqKfXYG3DbT0bYXcvWatsNFZKOI2EUk/gLL+onIahGZ4Uqb6k8SWhVu/IZ9LUfSPHkaO17tz9FDF74E3NoFk2n9+1h2BjSl7v0zqRQW4fmwbhY/ZBQ5ra8lMuENxy+bstbCFyCkmvft1fZ9EppdAj8/7rnPyYZv4fOhEBoJd82DuucthZZydQ9/AzAMWFyMZR8AklxsT53L5ke9Ef8lMf5VGuZswz6+N9tWn+e/Y/t82v4+hn2BTag1ehZh4V6yN1YKgUP+h4loQMbXt5OelmJ1nIprzx+wY4HjzNFgL/tL0eYHwz6EiPqOgdQTB9y3bmMcM4FOuwOi4xzH2JfVeQel5OpFzJOMMRc8BVRE6gKDAS85ONb3dLz8bvYP+x47ftT/fhhrZvz9Yump6+fC5Bux1WhB44fmEF61ugVJ3SgojH393yYo8wjbP7lTp16wgjGw4DmoXAs63WV1msKdHcQ9DVNudkxe5qr8PMcsnfOfgdbDHGfPektX1nmUVR/+G8BjgP5GelDTdj0IvO9XdgS3on3C4zD7qbMXS9/4+0xCp93I8ZB6cPMPSDn4cBZH/ba9SGwymo4Zv7LyuzetjlPx7FgAe/+A3o96zaGHharREq58H/YnwKxHXRvEzc6AyTdAwifQ40HHvDgBwW6L6kkXLPgiMl9ENhRyG1qcBkTkcuCIMSaxmMuPFJEEEUlISdE/00uqWo1oWj22ADrfA8ve5fD7g1kz7ysazb2dw341yb/pe6gUaXVMt+p84zjWB3WgzbqX2LNltdVxKg5jHBPshdeHuFutTnNhrYZAr0dg1eeQWMrzOE4egs8GwfZ5jpkuBzzrOPa/nBDjhsOVRGQR8KgxJqGQ514CbgbygGCgCjDdGHPThdYbHx9vEhL+tkpVTGm/fULlef8gUPLYJ3UIuWcO1WvVtzqWRxw9sAe/CT055hdJw8eXYgv0/kNMy72kGTDlRhj6rmPenPLAnu+4wtTORXDbjJJdhetIEkwa7pjiePhn0PwST6V0iYgkGmMKHTn2+FeTMeZJY0xdY0xD4DpgYXGKvXJdtZ53sOuKb1gZMZCgu2b5bLEHqF6nAfv7vkbj/F3YFjxrdRzfZ8+HX16AyKYQe53VaYrP5gdXfwQR9ZyDuAeL97qdv8LHl0J+Ltw+y2uL/YW4eljmVSKSDHQDZorIHOfjdURkljsCKte0iO9PpwcnUyO6kdVRPK5N32sdp7Ivf59TG/Tj51Ebv4MjmxyHPXrBjKolEhIBIyY5+uKnFmMQd83X8OXVUKWOYwK0Ou3LIqVHuKVLx1O0S0eVWG4WaW/2goxD2Ef+RvU6DaxO5Hvy8+DdzuAfDPf+Vq76sP9i0w+Ovfy4W2HIW39/3hjHzJ+LXoRGveHaLxxfFl7O0i4dpcpUQDAnL/+AEJPFwYm3Ys/PtzqR71n7FaTtcFyYvLwWe4BWQ6Hnw7Bq4t+nJ8nLgR9GO4p9u+vhxm/LRbG/kHL8v6VU4RrExLGu7ZO0zV7Niq+0P9+t8rIde73RHR3T/pZ3/f8JTS+GWf+Avcsdj2Wlw1fDHRdH7/OE43BOf++4NoSrtOArn9R52IOsqtSbjtvfYesqz18prMJInAjp+xyFsoyvx+oRZwZxw+s6+vOTE+GTy2D3b47Le/Z70jfep5MWfOWTxGajyZ2fkGarSp0FoyH7pNWRyr+c07DkVWjQAxr3szqN+4RUdZyJm50BH/WH9GS4cRp0uNHqZG6nBV/5rPBqUUTe8jmVT+93/MmuXLNiAmQchv5P+9ReL+C4BOHVH0K9rnDHbGjiQ19oBWjBVz7Nv1EP6PM4rP2a1TM+sDpO+ZV1An5/w9Hf3aCb1Wk8I2Yw3DkHara2OonHaMFXPs/e8xGSAlrTbOUz7N+pE7aWyrL3IPOYo+9elVta8JXPs/kHEHHTZ9jFxqmvbiE3xw2zJVYkp9Pgj3cg5nKo08HqNMoFWvBVhVC7QXO2d3mB5nlbSZj4mNVxypff34ScDN279wFa8FWFETfwdhKqDKBD8iTSUw9bHad8OHkYln8AbYc7phhW5ZoWfFWh1Bz0OMGSS1jSZKujlA9LXoP8HOj7hNVJlBtowVcVSr2YTtCgJ7aEjxwzPqqiHd/nmDe+w40Q2cTqNMoNtOCrCicz7k44vpcNi6ZaHcW7LX7F8W9vHfPwFVrwVYXj3+pyDhOJWT7B6ijeK3UHrJ4EHW93zB2vfIIWfFXhBAQEsrPhCNpmr2LPljVWx/FOi14Cv0DHJQGVz9CCryqk5gNHk2P8OTSvkHnQK7rDm2D9NOhyD4TVtDqNciMt+KpCiqxZl7URF9EmZSYn09OsjuNdfnkBgsKgxwNWJ1FupgVfVVjhfUZTSbIwa762Oor32L8KNs+AbqMhtJrVaZSbacFXFVbzuD4QHU+VdZ+C3W51HO/wywsQUg263md1EuUBrl7EfLiIbBQRu4gUeg1F53IRIjJNRDaLSJKI+Oh0e6rc6TwSUrexL/Fnq5NYb89S2D4fej4IwVWsTqM8wNU9/A3AMGDxBZZ7E5htjIkB2gE6ZaHyCrkxQ0glnNRf3rY6irWMgYXPQeWa0Oluq9MoD3Gp4BtjkowxW863jIhUAXoDHztfk2OMOe5Ku0q5S0BQCFvrXk3sqWUc2LXZ6jjW2fkL7Pkdej0KgaFWp1EeUhZ9+I2BFOBTEVktIh+JSKWiFhaRkSKSICIJKSkpZRBPVXSNB47BjrB3TgU9RNMYWPAchNeDjrdanUZ50AULvojMF5ENhdyGFrMNfyAOeN8Y0wE4BRQ5E5MxZoIxJt4YEx8VFVXMJpQqvZrRjVkX1ouWh74n81QFvPbtlllwYBX0eQz8g6xOozzoggXfGHOxMaZNIbcfitlGMpBsjFnuvD8NxxeAUl4juMcowjnF4d+/sDpK2bLbYeELUK0JtLvB6jTKwzzepWOMOQTsE5EWzocuAjZ5ul2lSqJll0vJj2pNwx1fObo4KoqN0+HIRuj3FPj5W51GeZirh2VeJSLJQDdgpojMcT5eR0RmFVh0DDBJRNYB7YEXXWlXKXcTmw2/rvfA4fVk7vjN6jhlIi83h5QZ4zgY1Ij/HWzDrPUHzz6XlatTR/siMV68NxMfH28SEhKsjqEqCJNzilMvNWdrpY7EPfqj1XE8bvk3r9Jl43M8Ynuc6Znt6NY4kq/u7gpAv1cXkXIym9rhwdSOCKFOeDCdG1VjWFxdAPamnqZ6WCChgfpXgbcRkURjTKHnRen/llJOEliJDTWHEn/wa47s30WN6EZWR/KY7NPpNN34NkkBrXn1ySd42cDp7D/36m/p1oA9qac5cDyTQyey2HTgBPl2w7C4uhhjuOSNX8nKtRMeEuD4UggPZnBsHa7p6Hh+6c5UaoeHUDs8mOAAPwvfqSpIC75SBdS7ZAy2iV+x4+e3qHHX61bH8ZigFe8TxHHSB32G2GwEAOGhf/bw3t7j7192Z3oD7AZeHhbLgfRMDqVnceB4FgfTM0k5mQ3A8dO53PDh8rOvqxoaQO3wEO7u3YirOtTldE4eczcedn5RhFAzPIggf/1SKAta8JUqILpxS9ZU6krz5G/JznqBoGDfOwnJnDyE/P4WtBpK4w79iv06EQHAzyZc2SG6yOVCg/z46q4uHEx3fBEcSM/i4PFMgp1FfffR0zw4Zc1fXlO3aghvXd+BuPpVS/6GVLFpwVfqHH5dRhL5y+0kzJ1I/JBRVsdxu4SJTxCXm4Wt/78QD6w/yN+P7k2rF/l80xqVWfBIHw4ez+JAeiYHj2exfFfq2S8E5Tla8JU6R+ueQ8lY2Yi4Q98AvlXw925dQ4eUH0iIuoou1ZtakiHQ30aTqMo0iapc4NFmZ3+y2w02mye+ipROj6zUOWx+flTuNQrbgUTYn2h1HLdK/eEpsgii6fB/Wx3lb4wx/PP79fzzhw1WR/FZWvCVKky768n1C2XD969ancRtNi+fQ4dTv7O+0W1E1qxrdZy/EREqBwXw1fK9fL1ir9VxfJIWfKUKE1yFxKoDaX5kLmmHk61O4zpjsM3/F0eoRvvh/2d1miL949IW9GpWnWd+2MiqvcesjuNztOArVYTaA8YSKHls/fk9q6O4LulHmudu5kTXRwmpFGZ1miL52YS3r+9ArfBg7v0ikSMnsqyO5FO04CtVhAYt2rM+KI5GuyeTl5tjdZxSM3k5MP9ZiIqh6YB7rI5zQRGhgXxwc0fsxrDtSIbVcXyKFnylziMv/i5qksq6BV9ZHaXUVnz7OqTtILPPv8rNBGkta1dhyWP96XGewztVyWnBV+o8YvuNIDWgFk13lc+Cn3HiGM2S3mFjYFuCWw20Ok6JhAT6YYzhy2V7mLpyn9VxfIIWfKXOw8/fn8i+91Hl8HI4vNHqOCW2YepzVOMEAZc9j9jK36+7MTBn4yH++f0GHcR1g/L3CVCqrHW4GeMXTPKcN61OUiIpB3YTu+9LEsP60Tyur9VxSsXmHMStGR7EqC8TOXJSB3FdoQVfqQsJrcbKKhdRbcf3pKeVn+ss753+NP7kUeuqF6yO4pKI0EAm3BzPicw87vtyFTl5dqsjlVta8JUqhsh+9xMq2SSVl0M0j2wmLnUGh1vcRHTj1lancVnL2lX4zzWxJO49xm/by8+Xbmn8ujXl7Myj7qYFX6liaBLbnaSA1tTbPgl7vvdfDcrMfwYJrEy9oc9YHcVthrSrw9wHe9M/pqbVUTzizPTTPZtWZ//xTI+0oQVfqWI61f5Oos1h1v/6jdVRzmvjH7OQrbNJ7TAaKkVaHcetmtV0nDS2bGcqa/YdtzaMG20/cpIRHyzjYHomfjahfb0Ij7SjBV+pYmo34CZSqEbE+olWRymSsdsJWPgMh4mkUu/7rY7jEXn5dp6cvt5xJq4PDOIu35nKsPf+YOfRDFIzPHuCn6sXMR8uIhtFxC4ihV5D0bncQ87lNojI1yIS7Eq7SlkhIDCIiN730ODYH3B0u9VxCrV69mc0z9vKnnYPERxa+cIvKIf8/Wy8e0McxzNzGD2pfA/i/rBmPzd/vILqYUF8d18P2kSHe7Q9V/fwNwDDgMVFLSAi0cBYIN4Y0wbwA65zsV2lLBHQ6XawBZC9dILVUf4mJzuLmitfZpetIR2v8K15/M/Vqk4V/nN1LCt3H+P5mZusjlMqP649wAOT19C+fgTTR3WnXjXPX13NpYJvjEkyxmwpxqL+QIiI+AOhwAFX2lXKMmE12Rp1MTmJn3Pq5HGr0/zFxh9fJ9oc5kSvp/HzLx9TKLhiaPto7u7ViM+X7mHJtvJ35E7fFlGM7teEL+7sTERoYJm06fE+fGPMfuBVYC9wEEg3xswtankRGSkiCSKSkJJS/v4Tle+zd7qbMDLZMOsDq6P8KSud9jsnkFqjK7F9hlmdpsw8flkM/7u2HT2alI85dzKy83jp5ySycvOpEhzAPy6NKdMLuF+w4IvIfGff+7m3ocVpQESqAkOBRkAdoJKI3FTU8saYCcaYeGNMfFRUVHHfh1JlpkVcP7b5N6PWli8wdu/oP7YveQPJTCPyypfL5RQKpeXvZ2NYXF1sNmFf2mmPHb/uDodPZHHt+KV8tGQXK3alWZLhgp8MY8zFxpg2hdx+KGYbFwO7jDEpxphcYDrQ3ZXQSllJbDaOt7mNBvZ9bPz9J6vjcDh5B7m/v8ORBkOgTger41giOy+fER8sZfRXq8jN944v4YI2HzrBle/+zp7UU3x8azy9m1uzM1sWuwJ7ga4iEioiAlwEJJVBu0p5TNtLb+cYVchbZn23zp5p/0SMnby+3nslK08L8vfj8YExrNiVxgszvau8LNuZyvD3l2I3hqn3dqNvixqWZXH1sMyrRCQZ6AbMFJE5zsfriMgsAGPMcmAasApY72zT+w5xUKoEgkMqcbrNjbQ7vRSOW3f91V0bl9Px2M+sqjWcOo1iLMvhDYa2j+auno347I/dfJPgPdMp16wSTJvocL67rwet63j2sMsLkTOn83qj+Ph4k5CQYHUMpQqXngxvtIXuY2HAs5ZEWPufATTK3IgZs5rwSN+ccqAk8vLt3PLJChL2HGP6qO4eP669KMYY5mw8zKWta+Lo2Cg7IpJojCn0vKiKM7qjlLuF1yW13gAyln5M1umyvxTfrhWzaJe5gk1N7tJi7+TvZ+OdG+K4oXN9GkdVsiRDTp6dR6au5d4vE5m36bAlGYqiBV8pFxxqcSuV7SdZN+eTsm3YbqfR6v+QXakO7a95vGzb9nLVKgUybkhrQgP9OZ2TV6aDuOmZudz26Qqmr97PIwOaM6CVd30Ra8FXygWtug1kl60B1TZ8VqaHaOatnwYH1xB0yTMEh1izJ+vtTufkMey9P8psEHf/8UyGj/+DlbvT+N+17RhzUbMy7865EC34SrlAbDaOxNxC0/wdbElYUCZtZmedJuX7/yOtcnNoe22ZtFkehQb606NpdT77YzfTEpM93t6ulFMczchh4u2dGRZX1+PtlYYWfKVc1Gbg3ZwglFNLyubiKKu/fZXa5gj7Oz8FFegkq9J4cmAM3RpH8tR361mXfNwjbZyZu75ns+oseawf3Zt671m/+mlRykWVwsLZU+8qOmQshpOHPNpW+rEUWm4bz7rgjrTtfZVH2/IFjkHcDkRVDuLeLxI5muHeM3G/WLaHvv/9hcVbHdPAVAry7jmMtOAr5QZtr3wEm8mHhE892k7S1GcJM6epNOh5j7bjSyIrB/HBzR2pFR5MtpumUrbbDS/NSuLp7zfQq1kUHRtUdct6PU0LvlLuENmEvCYXk7nsI3KyPXNRjozDu+hwYDKJEZfQJFZnJymJNtHhfDuqO9ERIbh67lFWbj5jJq/mg8U7ualrfSbc3NHr9+zP0IKvlJtsrHsdIdlHWTfvc4+sv/If/yHA30aDa170yPp9nYhwOiePe79MZPqq0g/iztl4iJnrDvLkwBieG9oGf7/yU0bLT1KlvFzb3lexT+pQea37u3VyktfA2snYutxDjXpN3b7+iiLAz0Z6Zi5PTl/Phv3pJXptnvN4/iHt6jBjTE/u6dPE6w67vBAt+Eq5ic3Pj/3NbiQmdxPb1/7m1nVv+fJhTvtVhl4Pu3W9FU2A80zcyEqB3PNFIqnFHMRdtfcYF/3vV5IOnkBELJuywVVa8JVyo1aDRnHaBJH2y7tuW+f6xd/RNiuRdY3vhpDyMTjozapXDuKDm+NJychm9Ferzu65F2X2hkNcP2EZAEH+5btklu/0SnmZKhGRrK8+kPbH52HPOOry+uz5+YT++m8OSA06XP0PNyRUAG3rhvPSVW3ZdjiDPWmni1zu4992MWpSIq3qVGH6qO40jirfF4bXgq+Um8Vd8xiB5GJb84XL61o1cwJN8ndyIO5RgoI9f5HriuTqjnVZ+GhfmhRRxH9Ys5/nZmzista1+PrurkRWDirjhO6nBV8pNwuo3Roa9SZ/+Ufk5eaUej0mN5N6q19ju18T4gbd5caE6ozwkADsdsMb87f+bRB3YJvaPH9lG969IY7ggLK77qwnacFXygP2NLkJv5PJrFs4tdTrkBUfUtOkEDL4RWx+vlFwvNGJrFymrNzHPV8ksvXwScZ8vZq0UzkE+tu4qWsDbLbydSTO+WjBV8oDortcxSGiCFz1Yalen3XiKGbJq9B0ANFxl7k5nSooIjSQD27uSEpGNpe+sZh5mw6x+dAJq2N5hBZ8pTzAPyCQXY1G0CZ7DbuTSn7VtjWT/onJOkFu/2c8kE6dK7ZuBP+9JpYWNcOYPLIb3Zt47wRornD1mrb/FZHNIrJORL4TkYgilrtMRLaIyHYRecKVNpUqL2IG3U+2CeDw/LdL9LoDu7cQd+gbEiIGElCnrYfSqXMNbR/N7Ad7075ehNVRPMbVPfx5QBtjTCywFXjy3AVExA94FxgItAKuF5FWLrarlNerGlWbdVUvpu3Rn8k8cazYrzsw/SnsCA2u0QnSlHu5VPCNMXONMXnOu8uAwmb97wxsN8bsNMbkAJOBoa60q1R5Ue+yBwmVbEI2TS7W8tvX/kb8ifmsjr6BmnWbeDidqmjc2Yd/B/BzIY9HA/sK3E92PlYoERkpIgkikpCSkuLGeEqVvVoxXaFuZ1gxAS50CURjMHOe5hhhtL5W++6V+12w4IvIfBHZUMhtaIFl/g/IAyYVtopCHityflJjzARjTLwxJj4qKqo470Epr5Yeezuk7WTd4unnX3D7ApqdXkVWt4epEhFZNuFUhXLBSZyNMRef73kRuRW4HLjIFD7RdDJQr8D9usCBkoRUqjwLaTeMo7OexCz/APpeU+gy+Xl5yNynsVVtSO2L7i/jhKqicPUoncuAx4EhxpiiJqRYCTQTkUYiEghcB/zoSrtKlSeBQcFsqzectqdXsn/HhkKXWTXjfWwpmzja5QnwDyzjhKqicLUP/x0gDJgnImtEZDyAiNQRkVkAzkHd+4E5QBIw1Riz0cV2lSpXmg68n3xs7Jv790M0s05nUH/N62z1b061TtdakE5VFC5dl8sYU+iVGIwxB4BBBe7PAma50pZS5VlUnYYkVOlDq8M/cjrjFUIr/zmf+upp/6EbqRzt97ZOoaA8Ss+0VaqMhPe5jyqcxrb+z/l1jqUcpM2OD1kT0pXWPQZbmE5VBFrwlSojzTpeDLViCV79CTiPb9j/43OEkkXEkBcsTqcqAi34SpUVEfI7jYQjm9i5cg6k7aL1/qmcaDmChi3jrU6nKgAt+EqVodyWV3GMMI7/+i65855FbP5UHTTO6liqgtCCr1QZCg6tzObaV9IuYwkBSd+xL+YOqFLb6liqgtCCr1QZa3jZGADSqELEgEctTqMqEpcOy1RKlVztBi1YFvMYITWb0i68mtVxVAWiBV8pC3S9/imrI6gKSLt0lFKqgtCCr5RSFYQWfKWUqiC04CulVAWhBV8ppSoILfhKKVVBaMFXSqkKQgu+UkpVEFL4ZWi9g4ikAHtK+fLqwFE3xinPdFv8lW6Pv9Lt8Sdf2BYNjDFRhT3h1QXfFSKSYIzROWfRbXEu3R5/pdvjT76+LbRLRymlKggt+EopVUH4csGfYHUAL6Lb4q90e/yVbo8/+fS28Nk+fKWUUn/ly3v4SimlCtCCr5RSFUS5L/ji8JuIDCzw2LUiMltELhORLSKyXUSesDJnWTjPtlggIr+ISJKIbBSRB6zMWVbO99lw/uwnIqtFZIZ1KcvGBX5PIkRkmohsdn5GulmZtSxcYHs85Pw92SAiX4tIsJVZ3coYU+5vQBsgCQgGKgHbgCbADqAxEAisBVpZndWibdEDiHM+HwZsrQjb4nyfDedzDwNfATOszmnltgAmAnc5lwkEIqzOauH26AHsAkKcy0wFbrM6q7tuPnGJQ2PMBhH5CXgcx3/c50ANYLsxZieAiEwGhgKbLAtaBgrbFsaY3ws8f1JEkoBofHxbQJHbY4eI1AUGAy/gKPw+r4jfkxSgN3Cbc5kcIMeqjGWpiO2xG8elX0NEJBcIBQ5YFtLNfKLgOz0LrMLxYY0HrgD2FXg+GehiQS4rnLstzhKRhkAHYHnZx7JMYdvjDeAxHH/xVCTnbouWOIr+pyLSDkgEHjDGnLIuYpn6y/YwxmSLyKvAXiATmGuMmWtlQHcq9334Zzg/oFOAL4wx2YAUtljZprJGIdsCABGpDHwLPGiMOWFVvrJ27vYQkcuBI8aYRIujlblCPhv+QBzwvjGmA3AK8PnxrjMK+WxUxdET0AioA1QSkZuszOhOPlPwnezOGzj26OsVeK4uPvSnWTEU3BaISACOYj/JGDPdslTWKbg9egBDRGQ3MBnoLyJfWhXMAuf+niQbY878xTcNxxdARVJwe1wM7DLGpBhjcoHpQHfLkrmZrxX8glYCzUSkkYgEAtcBP1qcyRIiIsDHQJIx5n9W57GaMeZJY0xdY0xDHJ+LhcYYn9mLKwljzCFgn4i0cD50ERVgbOc89gJdRSTU+XtzEY6BXZ/gS334f2GMyROR+4E5gB/wiTFmo8WxrNIDuBlYLyJrnI89ZYyZZV0k5UXGAJOcO0Y7gdstzmMZY8xyEZmGo18/D1iND023oFMrKKVUBeHLXTpKKaUK0IKvlFIVhBZ8pZSqILTgK6VUBaEFXymlKggt+EopVUFowVdKqQri/wFrNmEFs3c9aQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "outcomes = ['Y'+str(x) for x in periods]\n",
    "data.loc[treated_unit,outcomes].plot(linestyle = 'dashed', label = 'Untreated')\n",
    "data.loc[treated_unit + 'treat',outcomes].plot(linestyle = 'solid', label = 'Treated')\n",
    "plt.title(\"True Effect on Treated Unit\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# Untreated represents the actual counterfactual (not synth control)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8b57f1b",
   "metadata": {},
   "source": [
    "## Measure Treatment Effect with Synthetic Control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3a50b905",
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = np.matrix(pretreatment_data.loc[treated_unit + 'treat']).T\n",
    "X0 = np.matrix(pretreatment_data.loc[donor_units]).T\n",
    "X0_all_periods = np.matrix(data.loc[donor_units]).T\n",
    "\n",
    "outcome_cols = ['Y'+ str(x) for x in pretreatment_periods]\n",
    "Y_donors = np.matrix(pretreatment_data.loc[list(donor_units), outcome_cols]).T\n",
    "Y1 = np.matrix(pretreatment_data.loc[[treated_unit + 'treat'], outcome_cols]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f34411c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_features = len(features) * len(pretreatment_periods)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "55b6cf58",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find V such that W(V) minimizes the mean squared prediction error\n",
    "V_star, estimated_weights = find_best_V_w(X0, X1, Y_donors, Y1, n_features = n_features)\n",
    "\n",
    "synthetic_unit = X0_all_periods*np.matrix(estimated_weights).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a23696a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD5CAYAAAAk7Y4VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABHVElEQVR4nO3dd3zN1x/H8de5N5vESKJWiLZGZQhir9jU3qtGq6W1aZWqFi2dthalqNWBVtFSm1ilsWNvYkYikSHj3nt+fyT1QxMJd3wzzvPxuI/k3nu+3/O+N8kn3/sd5wgpJYqiKErOp9M6gKIoimIbquAriqLkEqrgK4qi5BKq4CuKouQSquAriqLkEqrgK4qi5BJ25iwshCgI/AJ4A5eBzlLKe2m0yw98D/gCEnhDSrkvo/V7eHhIb29vcyIqiqLkKgcPHrwrpfRM6zmzCj4wGtgqpfxCCDE69f6oNNrNAP6SUnYUQjgALplZube3NyEhIWZGVBRFyT2EEFfSe87cXTptgMWp3y8G2qbRuRtQF1gAIKVMklJGmdmvoiiK8ozMLfgvSClvAqR+LZRGmxeBcGCREOKwEOJ7IUQeM/tVFEVRnlGGBV8IsUUIEZrGrU0m+7ADKgFzpJQVgThSdv2k118/IUSIECIkPDw8k10oiqIoGclwH76UslF6zwkhbgshikgpbwohigB30mgWBoRJKfen3l/FUwq+lHIeMA8gMDBQDfSjKM8hOTmZsLAwEhIStI6iWImTkxPFixfH3t4+08uYe9B2LdAb+CL165onG0gpbwkhrgkhykopzwANgZNm9qsoylOEhYXh6uqKt7c3Qgit4ygWJqUkIiKCsLAwSpUqlenlzN2H/wXQWAhxDmiceh8hRFEhxPpH2g0GlgshjgEBwGdm9qsoylMkJCTg7u6uin0OJYTA3d39mT/BmbWFL6WMIGWL/cnHbwCvPnL/CBBoTl/P5NoBuLwLvOuAV1WbdasoWYkq9jnb8/x8zd2lk/Vc2gVL22EyGdHZOULvtaroK4qikAOHVoi8sJP+nvlZ5eoCxqSULX1FUWwqIiKCgIAAAgICKFy4MMWKFXt4PykpySJ9eHt7c/fu3Yf3d+zYQcuWLZ+6TFRUFLNnz7ZI/+asb/z48UyePNliOTIrxxX8/KUbcU+vZ7mbKya9XcpuHUVRbMrd3Z0jR45w5MgR3n77bYYPH/7wvoODAwaDQZNcTyvQRqPRouvLinJcwdeVrE51j3ZcdLDnr2LV1O4cRcki+vTpw4gRI6hfvz6jRo36z1aur68vly9fBmDZsmVUrVqVgIAA+vfv/8zFePz48bzxxhsEBQXx4osvMnPmTABGjx7NhQsXCAgIYOTIkezYsYP69evTvXt3/Pz8MBqNjBw5kipVquDv7893330HQGxsLA0bNqRSpUr4+fmxZs2aNNcH8PXXXz9cfty4cQ8zTZo0ibJly9KoUSPOnDnz3O+jOXLePnzgrRbjWLV8A2sTz/JqUhw4qAt7ldxrwroTnLxx36LrLF/UjXGtfJ55ubNnz7Jlyxb0ej3jx49Ps82pU6f45Zdf2LNnD/b29gwYMIDly5fTq1evZ+rr9OnTbN++nZiYGMqWLcs777zDF198QWhoKEeOHAFSdgMdOHCA0NBQSpUqxbx588iXLx///PMPiYmJ1KpViyZNmuDl5cXq1atxc3Pj7t27VK9endatW/9nfZs2beLcuXMcOHAAKSWtW7cmODiYPHny8PPPP3P48GEMBgOVKlWicuXKz/z+mStHFnxXR2cCnGqzSx/MkeBvCWj0vtaRFEUBOnXqhF6vf2qbrVu3cvDgQapUqQLAgwcPKFTov6O2pHWWyqOPtWjRAkdHRxwdHSlUqBC3b99Os7+qVas+PJd906ZNHDt2jFWrVgEQHR3NuXPnKF68OGPGjCE4OBidTsf169fTXN+mTZvYtGkTFStWBFI+GZw7d46YmBjatWuHi0vKuJGtW7d+6ntgLTmy4AMMaTSaveuD+fXsjwQ0HAnqFDUll3qeLXFryZPn/5+27ezsMJlMD+//e065lJLevXvz+eefP3Vd7u7u3Lt3Dw8PDwAiIyMffg/g6Oj48Hu9Xp/ucYNHM0kpmTVrFk2bNn2szQ8//EB4eDgHDx7E3t4eb2/vNM+Bl1LywQcf0L9//8cenz59epY4TTbH7cP/V7lCXpQ3lmSTczK3T2/WOo6iKE/w9vbm0KFDABw6dIhLly4B0LBhQ1atWsWdOykjtURGRnLlyn9H/A0KCmLp0qVAygHXZcuWUb9+/af26erqSkxMTLrPN23alDlz5pCcnAyk7IKKi4sjOjqaQoUKYW9vz/bt2x/meXJ9TZs2ZeHChcTGxgJw/fp17ty5Q926dVm9ejUPHjwgJiaGdevWZeo9srQcW/ABulYZSbxOx0+7v9I6iqIoT+jQoQORkZEEBAQwZ84cypQpA0D58uWZOHEiTZo0wd/fn8aNG3Pz5s3/LP/RRx9x/vx5KlSoQMWKFXn55Zd57bXXntqnu7s7tWrVwtfX9+FB1ke9+eablC9fnkqVKuHr60v//v0xGAz06NGDkJAQAgMDWb58OeXKlUtzfU2aNKF79+7UqFEDPz8/OnbsSExMDJUqVaJLly4EBATQoUMH6tTR5uxBIWXWHZ8sMDBQmjsBSsv51TGKaNZ03IxD/uIWSqYoWdupU6d45ZVXtI6hWFlaP2chxEEpZZojG+ToLXyA+iW6E2Zvx68bxmsdRVEURVM5vuC/06A/BQw6NkbvAWOy1nEURVE0k+MLvou9I9Wca3DQ2Y69O7LPFXGKoiiWluMLPsDAJuNwNEl+PbdE6yiKoiiayRUF39ujCBVMJdnpmMDVc2owNUVRcqdcUfABulUfRaJOx7Kdau4VRVFyp1xT8Bv51KVMogvbxBUSYtTk6IpibZMmTcLHxwd/f38CAgLYv39/xgul4ffff+fkyf/PihoUFMSznK59+fJlfvzxx4f3Q0JCGDJkSKaXT05OZvTo0ZQuXRpfX1+qVq3Khg0bMr38o558LZn1rK85Pbmm4AM09u7ObTs9y9ePy7ixoijPbd++ffzxxx8cOnSIY8eOsWXLFry8vJ5rXc9bJP/1ZMEPDAx8OHpmZnz00UfcvHmT0NBQQkNDWbdu3VOv1n2ap70WWwwZnasK/hsNBlIoWbA5KhgeGcNDURRSpgbdNSXlq5lu3ryJh4fHw/FsPDw8KFq0KFu3bqVdu3YP223evJn27dsDkDdvXj788EMqVKhA9erVuX37Nnv37mXt2rWMHDmSgIAALly4AMDKlSupWrUqZcqUYdeulONy6Q1tPHr0aHbt2kVAQADTpk17bKKU2NhYXn/9dfz8/PD39+fXX3997HXEx8czf/58Zs2a9fC1vPDCC3Tu3BmAn376CT8/P3x9fRk1atTD5TL7WoKCghgzZgz16tVjxowZbN26lYoVK+Ln58cbb7xBYmKi2T+LR5k1eJoQoiDwC+ANXAY6SynvpdFuOPAmIIHjwOtSymebfdcCHOzsqO1Snd+S97E5eB6Ng962dQRFsb0No+HW8ae3SbwPt0NBmkDo4AVfcHRLv31hP2j+RbpPN2nShE8++YQyZcrQqFEjunTpQr169WjQoAEDBw4kPDwcT09PFi1axOuvvw5AXFwc1atXZ9KkSbz//vvMnz+fsWPH0rp1a1q2bEnHjh0frt9gMHDgwAHWr1/PhAkT2LJlCwsWLEhzaOMvvviCyZMn88cffwApQyL/69NPPyVfvnwcP57y/ty793j5On/+PCVKlMDN7b/vxY0bNxg1ahQHDx6kQIECNGnShN9//522bds+02uJiopi586dJCQkULp0abZu3UqZMmXo1asXc+bMYdiwYen/HJ6RuVv4o4GtUsrSwNbU+48RQhQDhgCBUkpfQA90NbPf59a/+Sc4mySrzyzSKoKiZD0J0SnFHlK+JkSbtbq8efNy8OBB5s2bh6enJ126dOGHH35ACEHPnj1ZtmwZUVFR7Nu3j+bNmwPg4ODwcMu7cuXKDydDScu/nwoebbdp0yaWLFlCQEAA1apVIyIignPnzj0155YtWxg4cODD+wUKFMj0a/znn38ICgrC09MTOzs7evToQXBw8DO/li5dugBw5swZSpUq9XBMod69ez9cn6WYOzxyGyAo9fvFwA5gVBrt7ABnIUQy4ALcMLPf51a0QGGqGIuz1zGMsxf2U+alalpFURTbeMqW+EPXDsDi1inzQOsdoMP3Zs8Wp9frCQoKIigoCD8/PxYvXkyfPn14/fXXadWqFU5OTnTq1Ak7u5QyZG9v/3AI4acNZwz/H/r40XbpDW386Bb9k6SUTx22+OWXX+bq1avExMTg6ur6n2XT8yyv5d/hmW0xrpm5W/gvSClvAqR+/c8sBVLK68Bk4CpwE4iWUm5Kb4VCiH5CiBAhREh4uHXOpulRfRQGIViyY6JV1q8o2Y5XVei9Fhp8mPLVzGJ/5syZx7aujxw5QsmSJQEoWrQoRYsWZeLEifTp0yfDdWU0pPG/0hva+GnLN2nShG+++ebh/Sd36bi4uNC3b1+GDBnycPL1mzdvsmzZMqpVq8bOnTu5e/cuRqORn376iXr16j33aylXrhyXL1/m/PnzACxdujTD9T2rDAu+EGKLECI0jVubzHQghChAyieBUkBRII8QIt0xTKWU86SUgVLKQE9Pz8y+jmdS07c+/olO7OIiMbGRVulDUbIdr6pQ512LzAMdGxtL7969KV++PP7+/pw8efKxKQ179OiBl5cX5cuXz3BdXbt25euvv6ZixYoPD9qmJb2hjf39/bGzs6NChQpMmzbtsWXGjh3LvXv38PX1pUKFCmzfvv0/6504cSKenp6UL18eX19f2rZti6enJ0WKFOHzzz+nfv36VKhQgUqVKtGmzdPL4tNei5OTE4sWLaJTp074+fmh0+l4+23LHmc0a3hkIcQZIEhKeVMIUQTYIaUs+0SbTkAzKWXf1Pu9gOpSygEZrd8SwyOnZ9nGKXx56wfecg5iSOdZVulDUbSS1YdHHjRoEBUrVqRv375aR8nWbD088lqgd+r3vYE1abS5ClQXQriIlJ1aDYFTZvZrtu6NhuGVJNkWtROT0ah1HEXJNSpXrsyxY8cynKxEsTxzC/4XQGMhxDmgcep9hBBFhRDrAaSU+4FVwCFSTsnUAfPM7NdsOr2eoDzVuOAoWbN7odZxFCXXOHjwIMHBwY/NOavYhlkFX0oZIaVsKKUsnfo1MvXxG1LKVx9pN05KWU5K6Sul7CmltOzVBM+pb/NPcTWaWKtO0VQUJRfIVVfaPsm9QFFqGYty2OE+xy4e1DqOoiiKVeXqgg/wWs1RSOCHnZ9qHUVRFMWqcn3Br+DTiMBER/6W57kXa97VhYqiKFlZri/4AC1LdiJGL5irJjpXFIuIiIggICCAgIAAChcuTLFixR7e//cCpme1Y8cO9u7d+8zLeXt7c/fu3efqM6dRBR9o02gELyea2Bm1DZMaRVNRzObu7s6RI0c4cuQIb7/9NsOHD39438HB4bmGAn7egq/8nyr4gM7OgQYulbnuYGL5rsVax1GUHKlPnz6MGDGC+vXrM2rUKC5cuECzZs2oXLkyderU4fTp0wCsW7eOatWqUbFiRRo1asTt27e5fPkyc+fOZdq0aQQEBLBr1y7Cw8Pp0KEDVapUoUqVKuzZswdI+XTRpEkTKlasSP/+/W0yRk12Ye7gaTlGr+YTWbmmGevPLqBnvde1jqMoFvPlgS85HXnaoussV7Aco6qmNU7i0509e5YtW7ag1+tp2LAhc+fOpXTp0uzfv58BAwawbds2ateuzd9//40Qgu+//56vvvqKKVOm8Pbbb5M3b17ee+89ALp3787w4cOpXbs2V69epWnTppw6dYoJEyZQu3ZtPv74Y/7880/mzdP8sp8sQxX8VPncS1Df+AKr7cPZf+Ew1V6qqHUkRclxOnXqhF6vJzY2lr1799KpU6eHz/072UdYWBhdunTh5s2bJCUlUapUqTTXtWXLlsdmj7p//z4xMTEEBwfz22+/AdCiRYtnGvI4p1MF/xHdq7/H2oOj+CH4U6q99JvWcRTFIp5nS9xa/h0K2GQykT9/fo4cOfKfNoMHD2bEiBG0bt2aHTt2PDbo2qNMJhP79u3D2dn5P889bcjj3Eztw39EWd9XqZ2g5xBnuR0TpXUcRcmx3NzcKFWqFCtXrgRSxoI/evQoANHR0RQrVgyAxYv/f0ztyaGFnxza+N9/HnXr1mX58uUAbNiw4T9DHudmquA/SghalehAvE4w6y81Vr6iWNPy5ctZsGABFSpUwMfHhzVrUsZeHD9+PJ06daJOnTp4eHg8bN+qVStWr1798KDtzJkzCQkJwd/fn/LlyzN37lwAxo0bR3BwMJUqVWLTpk2UKFFCk9eXFZk1PLK1WXN45HQlP6DLokrc1Tmy8fWD2On1tu1fUSwgqw+PrFiGrYdHznnsnWnqHMAdeyPf71qqdRpFURSLUQU/DV2bfoKnwcDGcwu0jqIoimIxquCnwaVQaZoaPDjvEMX2c2oUTSV7ysq7axXzPc/PVxX8dHStMgQHk2Tx7s+1jqIoz8zJyYmIiAhV9HMoKSURERE4OTk903LqPPx0lAzoQIO/J7DN+QzXoiLwyu+udSRFybTixYsTFhZGeHi41lEUK3FycqJ48eLPtIwq+OnR6Whboi1/Raxh2uYvmdrpK60TKUqm2dvbp3uFqpJ7mbVLRwjRSQhxQghhEkKkeRpQartmQogzQojzQojR5vRpS7UavkelhCT+ub+ZhOTnG9JVURQlqzB3H34o0B4ITq+BEEIPfAs0B8oD3YQQ5c3s1zac89PcyZcoOwOzdy3XOo2iKIpZzJ3E/JSU8kwGzaoC56WUF6WUScDPQBtz+rWl9o3HUSzZwOaLaqJzRVGyN1ucpVMMuPbI/bDUx7IFh6L+tEh2I8z+Hn+e2q91HEVRlOeWYcEXQmwRQoSmccvsVnpaw9ale66YEKKfECJECBGSVc4w6BQ4EGeTiSV/T9Y6iqIoynPLsOBLKRtJKX3TuK3JZB9hgNcj94sDN57S3zwpZaCUMtDT0zOTXVhX4UrdaBZv5Kw4zbm7N7WOoyiK8lxssUvnH6C0EKKUEMIB6AqstUG/lqO3p23xlhgETNk2Tes0iqIoz8Xc0zLbCSHCgBrAn0KIjamPFxVCrAeQUhqAQcBG4BSwQkp5wrzYtlepwXvUjE/gaMxm4lJn5lEURclOzD1LZ7WUsriU0lFK+YKUsmnq4zeklK8+0m69lLKMlPIlKeUkc0NrwrUwLR3LEWtnYNquZVqnURRFeWZqLJ1n0Lz+KLyTktl+eXHGjRVFUbIYVfCfgZ13LdokOnPH/h4rju3SOo6iKMozUQX/WQhBu0r9yWsysezgTK3TKIqiPBNV8J+Re2BPWscmcUWcJvTWVa3jKIqiZJoq+M/KIQ9tizVBIvl65yyt0yiKomSaKvjP4ZV6IwiKf8DJuC1EPYjTOo6iKEqmqIL/PNxforV9KRL0Br4KVqNoKoqSPaiC/5wa1BnBy0lJ7Lq2FJPJpHUcRVGUDKmC/5x0ZZrS6YEdUfZRLDm8Tes4iqIoGVIF/3np9LTy70N+o5Gfjs3VOo2iKEqGVME3g2vVvrSLSeCGOMM/Yee1jqMoivJUquCbw6Ug7QrXRY9k8u45WqdRFEV5KlXwzVSqzjAaxcVzLn4b4bH3tY6jKIqSLlXwzVU0gPa6IiTrDXwevETrNIqiKOlSBd8CqtcYjE9iIntv/ozBaNQ6jqIoSppUwbcAnW87usRJ4uyimReyQes4iqIoaVIF3xLsHGnu0wMPg5EVJxdqnUZRspUkYxJrL6wlPjle6yg5nir4FuJU7U06xcQSoTtH8KVsN4Oj1Y399QPaza9H+JntWkdRsphlp5bx4e4PGbN7DCaprlq3JlXwLSVfcdp6VMNeSqbum691mixlzMqRrIn9g/MOkUzd8hZcO6B1JCWLiEmKYWHoQtyd3Nl6dSvzj6m/HWsydxLzTkKIE0IIkxAiMJ02XkKI7UKIU6lth5rTZ1ZWtOYgmsfGcTlxB2HRkVrHyRJGrJ/Juvi/aBgXz1tR0fyR15m5u2ZrHUvJIpaeXEp0YjTfNvyWFi+24Nsj3xIcFqx1rBzL3C38UKA98LSfkAF4V0r5ClAdGCiEKG9mv1lTqbp0kvkx6ox8FvyD1mk09/6G2WwOn09QXDxf3Y3mnXsxVEpIYH7iIf6+ekbreIrGohKiWHJyCY1KNMLHw4fxNcZTrmA5RgWP4lL0Ja3j5UhmFXwp5Skp5VP/cqWUN6WUh1K/jwFOAcXM6TfLEoKAagOolJDAgTurSDIYtE6kmQ83zmXDnTnUiX/AFwHv4vD6euwbjGWC3cs4SiPvbR5IXGKi1jEVDS0MXUh8cjwDAwYC4GTnxIz6M7DX2TN0+1Bik2I1Tpjz2HQfvhDCG6gI7H9Km35CiBAhREh4eLjNsllMha50jU0m0S6GWX+v0TqNJsZtms+6m99SM/4BXwWMJE/NAeBVFeq+i3f3nxgTI4m2u82bv3+kdVRFI3fi7/Dj6R9p8WILXi7w8sPHi+QtwpSgKVy9f1UdxLWCDAu+EGKLECI0jVubZ+lICJEX+BUYJqVMdwwCKeU8KWWglDLQ09PzWbrIGhxdaVimHS8YDPx2NvddeTtxywJW35hJ1YQEJlcaTd4a/R9v4Jyflm0X0PF+LKEJG/juwHptgiqamndsHkaTkQEVBvznuSqFqzCyyki2X9vOd8e+0yBdzpVhwZdSNpJS+qZxy/TmqxDCnpRiv1xK+Zs5gbMDh+pv0+1+DPf1F1mai8bK/3zrIlaETadyQgJTAj/CtdqbaTcsWYPhZXrxUlIS80M/5tzdm7YNqmjqeux1fj33K+1Kt8PLzSvNNt3Ldaf1S62ZfWQ226+qU3ktxeq7dIQQAlgAnJJSTrV2f1mCZ1na5ffnhWQTU44M48s1n2idyOomb1/Cz9emUjExkSmB48gX2Oup7d0afcSnxkJIkcBb6wapISlykTlH5qBDRz//fum2EULwUfWP8HH34YPdH3Ax+qINE+Zc5p6W2U4IEQbUAP4UQmxMfbyoEOLfz+q1gJ5AAyHEkdTbq2alzgYKvtyIFTduUDHxAcuiVtJraV8SkpO0jmUVM3YsY9mVr/FLTGRKlQkUDHwt44X0dvh1XsLwqDgidGcZ+ud0q+dUtHcx+iLrLq6jS7kuFM5T+KltneycmF5/Oo56R4ZuG0pMUoyNUuZc5p6ls1pKWVxK6SilfEFK2TT18RtSyldTv98tpRRSSn8pZUDqLefvuJUmCppMzLt1hx7R9zlsOkDQsh5cjLytdTKL+iZ4OYsuf8kriUlMrTYRj0rdMr9wwVJ0r/8Z9ePi2R25mDUn0z2Wr+QQ3x7+Fke9I319+2aqfeE8hZlSbwphMWF8sOsDdRDXTOpKW2spVRfsnLFHx+jIKMZHxBDPWdqu7syfZ0K0TmcR83b9xIKLX1A2KYlp1T6jUEDnZ16HLqAbHxashacxmYn73uV2bLQVkipZwamIU2y6some5Xvi7uye6eUCCwfyftX32Rm2kzlH1URD5lAF31q8qkLvtdBwLHRcRIc8Jfjx5g3ycp/Re/sxccePWic0y8LdvzD3/CReSkpmevWvKBzQ4bnX9ULbWUyI1ZGkj6LP6pEWTKlkJd8c+QZXB1d6+/R+5mW7lu1K25fbMvfoXLZe2WqFdFnHugvrGLt7rFU+zaiCb01eVaHOu+DbHvpuxrdiX9Zev0z5ZAO/XPmczis/zJYXZy3Z8wvfnvsU72QD06t/TRH/ZzpD97+c8lGr3QLeioohzLSPCduWWiaokmUcuXOE4LBg3vB9AzcHt2deXgjB2Opj8fPwY8zuMVyIumCFlNqSUjL36FzG7B7DjbgbJBgSLN6HKvi2YucIzb/Ao/Nylt67T8eYB5yKX0vQ0te4GpV9LjD7ce8KZpz9FC+Dkak1JlPcv5VlVlyiOv383qJiQgKrr0xVQy/kIFJKZh6eSUGngnQv1/251+Ood2Rq0FSc7ZwZun0o95NyzpSiyaZkPt77Md8e+ZZWL7biu0bf4WLvYvF+VMG3tXKv4vD2bsa5vMy4uxHEiRO0/rUjG88d1jpZhlbsW8m0M59Q1GBkWs0pePu1tOj6HYJG8wnFcZFJDNs8RA29kEP8ffNv/rn1D/38+5ldxArnKczUoKlcj7nO6ODRGE3Z/3TemKQYBmwZwO/nf+edCu8wqfYk7PX2VulLFXwt5CsOvf+gY8UBLL55GzciGLm7L18Fr9A6Wbp+/XsVk09PoJDBxNSa0yjlY4Uza/V2eHdexEf34oizu0HfXHD9Qk4npWTW4VkUzlOYTmU6WWSdlV6oxOiqo9l1fRffHvnWIuvUys3Ym/Ta0IuQWyFMrDWRAQEDSLl0yTpUwdeK3g4ajCWgywp+iUrEJymOpZc+pfuqcVluv/7a/av46tR43I0mptacTmmfZtbrrIA3zRt/Rcf7MZx4sJZ5/6gpI7OzHdd2cPzucd72fxsHvYPF1tu5bGc6lO7A/OPz2Xxls8XWa0snI07SfX13bsfdZk7jObR52cxjYZmgCr7WXgyiSP89LHIqR/uYWI7H/Ubjpa9lmfH0N+xfxWcnxpPfaGJKjRmU9Wli/U79OzPihfq8mJTM3OPj1dAL2ZRJmph1ZBYlXEvQ+uXWFl23EIIx1cbg7+nPh7s/5Ny9cxZdv7XtvLaTPn/1wV5nz5LmS6hepLpN+lUFPyvIWwinnqsZ5z+YDyOiiBahtF/Vhu0Xj2saa+P+lXx6YhyuJsnXNWZS3qexzfp2bT2dz+LsQMTRd90wNfRCNvTXpb84d+8cAwIGYK+z/D5pB70D04Kmkcc+D0O3DyU6MXtcw/Hz6Z8Zsn0I3m7eLH91+WOjhVqbKvhZhU6Hru4Iunb6lblRJpxFBO/u7Mm0PdqMNbdt/0o+PTEeZwlf1ZiJv08j2wZwyodPh4W8dy+ae7qTDN0wy7b9K2YxmAzMPjqb0gVK07xUc6v1U8ilENOCpnEz7iajdo3K0gdxTdLElJApTNo/iTrF6vBDsx/wdLHtiMCq4Gc1XlWp3n8Pi+3LUyY5noXnx9F31Yc23cIN/mcl406Mx17CF9VmUtGnoc36foxXVbpWHEBQXDy77i5k3Sk1F252sfbCWq7cv8KggEHohHXLTEChAMZUG8Oe63uYdThrbhgkGBJ4b+d7/HDiB7qV68aM+jOsctplRlTBz4qcC+D92krm+wyjZewDDsStpfXi9tyKuWf1rvcc+IWxx8ejk4LPq8+iiq9GxT6Vru5Ixtl742E0MH7vSMJjc8651zlVkjGJuUfn4ufhR32v+jbps1OZTnQs05EFoQvYeHmjTfrMrMiESPpu6suWK1t4v8r7fFD1A/Q6vSZZVMHPqoTAtdYAJrVbxZD7cEN3ga6/NGWXFffr/33gZz48PgGk4LNqM6nu08BqfWWa3g6Pjgv57F4cyfpIev0+WutESgZWnl3JzbibDKo4yKqnGD7pg6ofEOAZwEd7PuJMZNa4cO9S9CV6/NmDM5FnmBo0lZ7le9r0PXmSKvhZnK5oAG/128eXogxGXSzv7ujO3B2Wn0kr5MBPjDn+CQah45OqM6jlmwWK/b8KlKR60yn0i4omzLiLT7Yv0zqRko745HjmHZtHlcJVqFGkhk37dtA7MDVoKq72rlniIO7B2wd5bf1rxBviWdh0IY1K2vg4WBpUwc8OHPPStPdvzHmxP16GZGZf/or3fnwHk8kygysd3v8jo45/SoLQMaHKTIL8tN2Nkya/jvQr3pQKCYn8enkq+69mr9PwcosfT/9IZEIkQyoO0WRL1tPFk6n1p3In/g7vB7+v2UHc9RfX89amtyjoVJBlry7D39NfkxxPUgU/G/ENGsqCV3+h3gMdG5N302tBfcLv3zVrnccPLOP90InECz3jKs+goV8W2rJ/gkPLKXyW4IizTGLw5mHEJ6uhF7KS+0n3WRS6iDrF6hBQKECzHBU8K/BhtQ/Ze2MvMw7PsGnfUkq+P/49o3aNwt/Tn2WvLsPLNe1pHLWgCn42k794BWb0/ZvXDCU4bh9Bn18a8M+J55s399Tfixl5bBIxQs/YytNoWiELbtk/ysmNEu0XMj4ikgd2V3nj94laJ1IeseTEEu4n3WdwxcFaR6FDmQ50KduFRaGL2HDJNldrJ5uSmbBvAjMOzeDVUq8yr/E88jnms0nfmaUKfjakc3BhVN8/Ge3RmXt6I0P3D+aXP8Y/0zrO/v0DI0K/IFJvxwcVp9Cigvb7FzPFqwrNqg6jw/1YTsT9zvx/stYZGblVZEIkS08upXHJxrzi/orWcQAYVWUUlQpV4uM9H1v9IG5sUiyDtg7i13O/0s+/H1/U+cKiQ0lYirlz2nYSQpwQQpiEEIEZtNULIQ4LIf4wp0/l/7q1/Jjp1b+lgFHPZ3dX8eWC5pgS4zJc7sK+BQwP/ZK7ejveD5hMm4o2GC7Bkuq8y3vOL+OdbODb4+O4EHFL60S53oLjC0gwJjAoYJDWUR6y19szJWgKbo5uDN0+lKiEKKv0cyvuFr3+6sX+m/uZUHMCgysO1vRMnKcxdws/FGgPBGei7VDglJn9KU+oWj6IH7psxzfZg2V2YQxdWIOoq4fSbX9573yGh37NHb0dI/2/pmPFpjZMayE6PXk7LuDLe/HoRBx91g232AFs5dndjrvNz6d/puWLLXkx/4tax3mMh7MH04OmEx4fznvB72EwWXZgwtORp+nxZw9uxN5gdsPZtC/d3qLrtzRzJzE/JaXM8LOSEKI40AL43pz+lLR5urmztO82mjjUYaezib5/defk9mkg5WPtru2ew9ATU7huZ89wv6/oXNmKo15aW/4SlG8+hfciI4kSoQxZnzWvsMwN5h2bhwkT71R4R+soafLz9GNs9bHsv7mf6QenW2y9u8J20XtDb4QQLGm+hJrFalps3dZiq33404H3AbUZZiU6nY4p3WYz8MXRXLOz5+1L89m0uB0kpFyZen33Nww9NZ1rdvYM8fmS7oHWG9/EZvw60rVUS+rFPSA4fCHrTv2jdaJc51rMNX479xsdSneguGtxreOkq13pdnQr143FJxfz58U/zV7fyrMrGbxtMCXdSvJjix8pU6CMBVJaX4YFXwixRQgRmsYtU4M3CyFaAneklAcz2b6fECJECBESHp59pv7LKvrXfY1p9X9GmPLyPudZ8F1V7i5oyrCTs7hs58CA8p/Ru6oVJi/RiK7FZMYnu1DQZODjve+roRdsbM6ROeh1evr599M6SoZGVhlJ5RcqM27vOE5FPN/eZZM0Me3gND7Z9wk1i9bkh2Y/UMilkIWTWk+GBV9K2UhK6ZvGbU0m+6gFtBZCXAZ+BhoIIdK9VFJKOU9KGSilDPT0tO1IcjlFrVI+/NJlIy8aijI9vz3txDUuONgztOjrvFnNQnPQZhWOrnh0WMiXdyMw6u/S6/cPtE6Ua1yIusAfF/+gW7lu2aLo2evsmVJvCvkd8zNs+zDuJTzb2FSJxkTeD36fhaEL6VymMzMbzNRkADRzWH2XjpTyAyllcSmlN9AV2CalfM3a/eZ2hV0LsKJ0E965F40JmHY7gt4uMsPlsqXigVSr8S5vRUUTZgxm4o4ftU6UK3x75Ftc7F14w/cNraNkmruzOzPqz+Dug7u8tzPzB3HvJdzjrU1vsfHyRt6t/C5jq4/FTmdn5bSWZ+5pme2EEGFADeBPIcTG1MeLCiHWWyKg8vzsXqzLgNgkdl27Rb1kCd51tI5kPbVH8LZbefwSkllxaQr/hJ3XOlGOdiLiBJuvbKZn+Z4UcCqgdZxn4uPhw7ia4zhw6wBTQqZk2P7q/av03NCTE3dPMLneZPr49smyp11mREiZdbf6AgMDZUhIiNYxsrdrB+DyrpRi71VV6zTWFXWNa/Nq09HDDRPe7Oy5Chd7R61T5Uhvb3mb4+HH+avDX7g6uGod57l8eeBLlp1axme1P6PVS2nv6jxy5wiDtw1GIJjZYKamQ0ZklhDioJQyzeui1JW2OZ1XVajzbs4v9gD5vfB6dRqf3A0nQX+ZN9dM0jpRjnTo9iH2XN9DX7++2bbYA4wIHEGVwlWYsG8CJyJO/Of5jZc30ndjX/I55mPZq8uyRbHPiCr4Ss7i256mZdqnTAgfu5qFIZu0TpSjSCmZeXgmHs4edCvXTes4ZrHX2TO53mQKOhVk2PZhRDyIAFJe46LQRby38z18PHxY1nwZJdxKaJzWMlTBV3Ke5l8y0uhGiWQTM46N52Lkba0T5Rj7buzj4O2DvOX3Fs52zlrHMVtBp4JMrz+dewn3eG/neyQYEpj490SmHpxKM+9mzG8yn/xO+bWOaTGq4Cs5j6MreTssZEr4XYQuht5r31VDL1jAv1v3RfIUoWOZjlrHsZjy7uUZX3M8IbdDaLG6BSvOrqCvb1++rPsljvqcdQxIFXwlZypemXJ13mdkZCRR4ihDN3yjdaJsb9vVbZyIOME7Fd7JkiNBmqPliy3pXb43EQ8i+LjGxwyrPMzqk69rQZ2lo+RcJiOmxS0ZmHCB3U4ufFHze1qUfeqgrko6jCYjHdd1xGAysLrN6mx5DnpGpJTEJMfg5uCmdRSzqLN0lNxJp0fXfj4T7ydRwGTkw92jiIiP0TpVtrTh8gbOR51nYMDAHFnsAYQQ2b7YZ0QVfCVny1cc95bT+Tr8Nkb9HV5bPUbrRNlOsimZ2UdmU7ZAWZp4Z7O5E5THqIKv5Hw+ban2SifejL5PmGEHk9TQC89kzfk1XIu5xqCKg3Lkfu3cRP30lNyh2Ze8IwpSPtHA2guTib+4W+tE2UKiMZG5R+fi7+lPveL1tI6jmEkVfCV3cMyLQ+3hDLwXSbxdMlt/654y7ITyVCvOrOB2/G2GVBySbcePUf5PFXwl94gLp9aDREomJ/Ozq1PKGENKuuKT4/n++PdUK1yNakWqaR1HsQBV8JXcw7sOejtHutyP5ZiTIxuFh9aJsrTlp5YTmRDJ4EqDtY6iWIgq+Eru4VUVeq+jlZMXziYT0y/t0TpRlhWdGM2i0EUEFQ+igmcFreMoFqIKvpK7eFUlf+tZtI6N44ZhPxcibmmdKEtafGIxMckxDKo4SOsoigWpgq/kPsUD6WBfHJPOxGe7FmmdJsuJeBDBslPLaObdjLIFy2odR7EgVfCVXOmVWsOo/uABRyN+JyE5Ses4Wcr3x78n0ZjIgIABWkdRLEwVfCV3Kt+Gjgl6Eu3imb5vtdZpsoxbcbdYcWYFrV9qTal8pbSOo1iYKvhK7qS3p0HAGxRLNvDnuSVap8ky5h6diwkT71R4R+soihWYO4l5JyHECSGESQiR7jCEQoj8QohVQojTQohTQoga5vSrKJZgX6UvnWPjibK7yl9nD2kdR3NX71/l9/O/06lMJ4rmLap1HMUKzN3CDwXaA8EZtJsB/CWlLAdUAE6Z2a+imC+PBy2LN8bRJJl9YJ7WaTQ3++hs7HX29PPvp3UUxUrMKvhSylNSyjNPayOEcAPqAgtSl0mSUkaZ06+iWEqh2sNoGRvH1eR9XI0K1zqOZs7dO8f6i+vp9ko3PJzVBWk5lS324b8IhAOLhBCHhRDfCyHypNdYCNFPCBEihAgJD8+9f4CKjRT2pZ3Tixh1Jj7bmXtP0fzm8Dfksc/DGz5vaB1FsaIMC74QYosQIjSNW5tM9mEHVALmSCkrAnHA6PQaSynnSSkDpZSBnp6emexCUZ5fhdrDCXyQwOG7v5FkMGgdx+ZC74ay7do2evn0ylETdiv/lWHBl1I2klL6pnFbk8k+woAwKeX+1PurSPkHoChZQ9nmdEh0IN4ujm/+Xqt1GpubdXgW+R3z0/OVnlpHUazM6rt0pJS3gGtCiH8v2WsInLR2v4qSaTo9jSv3o7DBwLqzC7VOY1Mhx5ay98Ze+no1Jq9DXq3jKFZm7mmZ7YQQYUAN4E8hxMbUx4sKIdY/0nQwsFwIcQwIAD4zp19FsTTHyr3pGJPAXf0Vtl04pnUcmzBd+ZuZ+ybiaTDSadtMrm6ZS8S105ji74GUWsdTrEDILPyDDQwMlCEhIVrHUHKJm78NpGX0Tko61OW3HrO1jmN1vy7uwHjOMvZuJF1iYh97zoCeWJ0rCXb5SHbMj3QqiC6vO46uHrjk98Qlnycijwc4FwSXgilfnQuAPmdOcJ6dCCEOSinTvC5K/XQUJVWRusNo9vN61rOHG/cjKepWUOtIVpOQlMCCpDOUxEC72HhMOgfO+Q4j3ORK4v27GOMiEA8isUu8h3N0NPmiz1PgzhHciMVRJKe73mR7N0zOBdDncUefxx3hUhBc3FP/MRRI/er+/38SLgXB3jll9rHLu8C7Tsow1opVqIKvKP/yKE1bp3Ks1YXx+c5FzGr1rtaJrGb66pFccxCMsq+OQwM/8K5DWa+qpDU2pskkiYhL4mb0A45GPSA8IpLoyNvERoWTdD+c5NgI9A8icZMx5DfEUiAhhgL3YnHXXcBdd5R8xOIi49MPo3cAYzIIAXpH6L1WFX0rUQVfUR5RJehdKmwdQMidlRiMw7DT67WOZHHh0eGsj93OK0ZB9x7zwc7+qe11OoGnqyOero74F88PFAF8HmtjNEkiYhO5EZ3AzagHXIhOYHf0A25EJ3ArOoHwe/dJjI0kn7xPAWLJL2IoIGLx1MfRWHcUP04ipARjUsqWvir4VqEKvqI86qWGtNvgzHinOL77Zz0Dq7fSOpHFTf59IPfsBCO9+6PLoNhnll4nKOTmRCE3JwK88qfZxmiShMckciP6ATejErgZ/YCb0QlsvBVC2WsjcBDJIE2I4qrYW4sq+IryKCF4tco7fHNiJmtOzstxBf9M2Em2m05SLcmJVg1tO5uVXiconM+JwvmcoMSjz5Tn8J4XuPTXTNrrdyMPL0F4107ZxaNYlBoeWVGe4FypJ+1jk7ipv8yeKzlrnL/pfw0hWcBb1cZrHeUxFWs15U6jmUxO7oQ4tgL2ztI6Uo6kCr6iPMkhD21KtcJOSmbu+VbrNBaz69gm9tndIijJnWoBLbWO8x/9677I5fLv8KexGnLLODi3RetIOY4q+IqShhJ1h9Mk7gHnE3cRHntf6zgW8d3fH+MkJUOaztA6SpqEEHzVqQLzCr7HaVkC08rX4e45rWPlKKrgK0pa8pegVR4fknQmvtiZ/YdbWLHje446xtFMlqZUyQCt46TLxcGOb3rXZoR4n/vJAtOPXSEhWutYOYYq+IqSjpr13scnMZEDN3/BZDJpHee5mYxGfjz3DZ4GE0PbztE6Toa8CrrwYfem9E8cioy8hFzVF0xGrWPlCKrgK0o6dN61aJPkSpR9LAtCNmkd57nN+WM8FxyMtHGpSYGCRbSOkym1S3vQqFk7xib3QZzfDFvGax0pR1AFX1HSIwStqw2moNHI6tDsObZOzIMYfg//He8kyTvts+a++/S8WacU8X49WWpsDHtnwtFftI6U7amCryhPkadiN1rHGQnTXeKfsPNax3lmU38bxi176FqkEw6OLlrHeSZCCL5o788K94EcwAfT2sFw/aDWsbI1VfAV5WnsHGlbqi06JDODp2ud5pmEhV9jY+J+/BN0dGs+Vus4z8XZQc/sXtUYJUZw25QP00/dIeaW1rGs69oB2DUl5auFqYKvKBl4qd4IGsQlcDphN/fiYzNeIIuYum4QsTp43Wc4umw8JpBXQRcm9Qiib9IIkuKikD/3gOQErWNZx7UDyMUtkVs/hcWtLV70VcFXlIy4vkDLvBVI0Bv5akf2OEXz6IUQdooL1ErMS6PafbSOY7aaL3vQvllThia+jbgeAn8Mz5GTtJgOLwNDIgKJ/HcgOQtSBV9RMiGo4RjKJiax78ZP2eIUzW+2pQzt/E6dnDO5XN/apXCp0I5phg5w9Ef4O3seSE9P4v6FmA4tRUowoUsZNtq7jkX7UAVfUTJBV7wyrYwFiLCPZenhrVrHeapNIavZbx9BI0MR/Ms30DqOxQgh+Ly9H9sK9WaTrIbcNBbOZ+2fRaYYDcT/PhzHDcPZbfJjS+AcdA3HIqwwL4C5c9p2EkKcEEKYhBBpTqmV2m54artQIcRPQggnc/pVFC20qz6UfEYjvx3N2gN7LTw0CVeTZEjLb7SOYnFO9nrm9qrKBP0gLooSyJWvQ8QFrWM9v/hI4ha2weXIQhaZWmDs8jNNWnWHOu9aZU4Ac7fwQ4H2QHB6DYQQxYAhQKCU0hfQA13N7FdRbM7NvyMt4yWXdZc4cvOy1nHStHjTdE44JtJC70exImnNX5X9FcvvzOQetXgjYTixyRL5UzdIyIbjHYWf4cGcIOyv/80E/UAC+82hoU9Rq3ZpVsGXUp6SUp7JRFM7wFkIYQe4ADfM6VdRNKG3p92LHQHJzB1TtE7zH0nJyay8spAiySaGdMhZ+7efVOMld/q0qMdbCUMwRZyHX9/MXsMvnN1E8ncNiLt/j/fzTOLNwR/hVzyf1bu1+j58KeV1YDJwFbgJREsps+916kquVjZoOHXjEznxYBfRCU+Zp1UDM38fxRUHSft8jcibN+dOwP6vPjW9KRrQmHFJveDcRtj2qdaRMiYlcs9M5I+dOZvkzidFvuWTwX0plt/ZJt1nWPCFEFtS970/eWuTmQ6EEAWANkApoCiQRwjx2lPa9xNChAghQsLDwzP7OhTFNlwK8qprJeL1RqbuWKB1moeiYu/xx/1NlE6Efu0max3HJoQQfNbOj2NFOrJCNoLd0+D4Kq1jpS85AePqdxCbP2K9sQrLfeYx+c2WuDlZZprJzMiw4EspG0kpfdO4rclkH42AS1LKcCllMvAbUPMp/c2TUgZKKQM9PT0z2YWi2E7Txh/xUlISe8J+zDKnaH7920Ai7AQ9Sva22Dy12YGTvZ65r1Vmql1fjunKI9cMhBuHtY71XzG3MSxqgf7YT0xL7sCFet8wqXM1HOxse6KkLXq7ClQXQrgIIQTQEMhZ88YpuYqusC+vGj24bR/LyuM7tY7DxRvn2GI8RuUEezo0eVfrODZXNL8zM3pU482EIURKt5QrcWNuax3r/24cwfBdEMk3jjPIMAzvjp8wpFFZhAZz9pp7WmY7IUQYUAP4UwixMfXxokKI9QBSyv3AKuAQcDy1z3lmpVYUjXWsMQJXo4mVB6dpHYVp6weTIOCtyh/m2om/q73ozqBWNegZPwxDbASs6AmGRK1jQehvmBY05W5sEr34lB6vD6FdxeKaxTH3LJ3VUsriUkpHKeULUsqmqY/fkFK++ki7cVLKcqm7gnpKKbPAT0JRnl9Bv3Y0ewDndZc4dSdMsxz7T+5kl10Y9ZIKUCuwg2Y5soKe1UviU6k2QxP6wbX98OcI7YZfMJlg2yRY9TpHDCV4y+lrPh/QnRovuWuTJ5W60lZRnodOR5uXu2BCMnPrl5rFmL37AxykZFCjrHeaqK0JIfi0rS/XizVjrmwPh5fB/u9sHyQxNuUTRvBX/GIMYpLHVywc2IKXC7naPssTVMFXlOdUod4waj5I5mj8LuISbf+hdc3upRxyjKGJyZsyL1r+qszsKOUgbiUW2HVjt74acuMYuLjDdgGiriIXNsF0ej2fJPdky8tjWda/Lp6ujrbL8BSq4CvK83Jyo5lrZWLsjMzcadtTNE0mE0tOTqWgwcTQNnNt2ndWVySfM7N7BjLwQX+u23khV/SGyIvW7/jKPuS8+jwIv0yfpJFQfQBzewbi7JB1hqZWBV9RzNCy6QRKJCez88oym/b7/frPOOtooJVTFTw9vGzad3ZQxbsg77UOpFvsUBIMJvipOyTGWK/DQ0uQi1txI8GBVgkTqN+iGx+3Ko9el7UOoquCryhmsPMsTXPjC1x3iGGNjU7RjE+M59dbv+CVZGJQ+6w9kJuWXqtWglqBgbwRPxjT3bPwW/+Ug6mWZDTAhtGwdjD/UJ52yZ8wqkcrXq9VyrL9WIgq+Ipips613sfFZOKXf2xz4HT66hHcsIdOnm1wctb+QGBWJYRgQhsfEorX4nNjTzjzJ+yw4PwAD+7B8o6wfw5LeZXB4kPm92tME5/CluvDwlTBVxQzFfJ5lcYP9JzSXeJcuHXHBbx97xYb4ndTPlHQu9UnVu0rJ3C0S7kS93eHlvxh1wiCv4bQ38xf8d1zML8hpku7+MDQj8Vub7NqYB0qeOU3f91WpAq+ophLCNqW7o5BwDdbPrdqV1N+H0CUXtC77EB0ejur9pVTvODmxNyelRn1oBdnHX2Qvw+Am0eff4XntyDnNyA+JpLOCWO4VKI9v75dE6+CLpYLbSWq4CuKBQTWG0KVB8kcig8mPtk6p2ievHKM7ZyleoIzr9brb5U+cqrKJQvyYeuKdI8eRIzOLeUgbuwzDs4oJez7Frm8EzdFIRrHjKdEQEOWvFGNfC7ZY/wiVfAVxRIcXGjqWo0oOxNztlvnFM2Zm4ZhEPB2zWwwDHAW1L1aCRpX9aNbzBCMcXdTh19IytzChkRYMwg2jiHEsQaNosbQsWFNpnSuYPMB0MyRfZIqShbXrvkEiiUb2HFlqcXXvfPIevbZ36FBciEq+zW1+Ppzi/Gty+NUohIjk/vD1X2wYWTGwy/E3oHFreDIMpY6dqX7/QF80qkawxuX0WQANHOogq8oFuJQ0JvGpiJcdohlw4ndFl33vAPjcTFJhjSbadH15jaOdnrm9KjEHqe6LLXrAAd/gH++T3+Bm0dhXn1MN44yWv8uXyW254c3qtOxsnYDoJlDFXxFsaBudUbjbDLx89+WG1/np22zOeb4gOa6cpT08rXYenOrQm5OzHmtMhMfdOCwUzXkhlFwKY1puU+ugYXNSEg20ClpHLvsa/PrOzWp9bKH7UNbiCr4imJBRcs1JijBgVDdJS5H3DJ7fQaDgZ8vfEchg4kh7edYIKECUKlEASa08aNn1FtEOHrBit5w73LKkyYTbP8cVvQi3OVl6kZ9TFIhf1YPqEmZF7L3dQ+q4CuKJQlBm9LdSdIJZm8y/+Dq7LVjuehgom3eeuR3K2SBgMq/ulYtQdvq5egQPYQkgxGWtIWtn6Tsr9/5Bcfcm1P79gj8y5Xhl/7VKeTmpHVks6mCrygWVitoKBUTDByI302SwfDc64mJv8+ae39QKknyTnvtJ1rJiT5u6UOhkq/wdUIb5L1LsGsK8sputrm1pfX11+haozTf9QzExSFnXPOgCr6iWJqdA43dahBhZ2L+9ucfj33Kr4O5YyfoXqw7dvZZY3jdnMbBTse3PSpRwMGIJOWMGxOCkEhHPmrpw/jWPlluADRzqIKvKFbQqdmnvGAwsPXy842iefX2ZTYmH6RCgh1dm39g4XTKowq5OtGoeQcSpT0GqSNJ2lO3cVv61i6V7U67zIgq+IpiBU4FitHQVJxzDrHsOLXnmZef+sdA4nTwRoWRuXaeWlsqE9iQA3V/YIlTD662+onq9V7NeKFsyNxJzL8WQpwWQhwTQqwWQuRPp10zIcQZIcR5IcRoc/pUlOyie90xOJgkP+57tvF1Dp7ZR7D+CrUT3WhQvbuV0ilPqtewBW988A1lAxtpHcVqzN3C3wz4Sin9gbPAfz57CiH0wLdAc6A80E0IUd7MfhUlyytZtj61Ex05LC5zI/puppebvXMkOgkDg762YjolNzKr4EspN0kp/z0N4W8grcvPqgLnpZQXpZRJwM9AG3P6VZTsok3p10jQCeZsGJep9hv2r+CAYzSNjMXxKVvLyumU3MaS+/DfADak8Xgx4Noj98NSH0uTEKKfECJECBESHv6Mo9kpShbTIGgwPgkm9sXuwmA0PrWtyWRi0dEvyW80MbT1tzZKqOQmGRZ8IcQWIURoGrc2j7T5EDAAy9NaRRqPpTtakZRynpQyUEoZ6OnpmZnXoChZl96Ohm41uW0vWbzt6ZONL9k4hVOOSbSwD6BIoZdsFFDJTTIs+FLKRlJK3zRuawCEEL2BlkAPKdMcdi4MeHSW5eKAdacFUpQspEeLSXgYjGy6lP4omglJCawIW0rRZBOD28+2YTolNzH3LJ1mwCigtZQyPp1m/wClhRClhBAOQFdgrTn9Kkp24uJWiHrSi5OOcew5vTfNNrN+H8k1B0mHAs3JkyefjRMquYW5+/C/AVyBzUKII0KIuQBCiKJCiPUAqQd1BwEbgVPACinlCTP7VZRspXvdD7GTkh/3TPrPcxH37/JHzHbKJsKbbb/QIJ2SW5g1QISU8uV0Hr8BvPrI/fXAenP6UpTsrEyZutTY7kyIwxVuR4fzQr7/H5+asnogkXaC4d5vqXlqFatSV9oqio20erkn8TrB/PX/P0Xz7LVTbDGdIDDBkbYNh2iYTskNVMFXFBtpXn8QZRIlu+N2Y0o9RXPGX0NIEtCv6scap1NyA1XwFcVWdDrqu9Xmur1k+bbZ7Dm+hT32N6mX5E6Niq21TqfkAmqHoaLYUO8Wn7JiZRAbLy1Dd0ngaC8Z1FiNda/Yhir4imJDrq6e1JYl+cPhKlIIOicVpnSpSlrHUnIJtUtHUWysT/kO6AFPg4ERd0Lh2gGtIym5hNrCVxQbK2OIYtzdSIoZDOQxGODyLvCqqnUsJRdQW/iKYmvedWibYKJKogH0DuBdR+tESi6htvAVxda8qkLvtSlb9t511Na9YjOq4CuKFryqqkKv2JzapaMoipJLqIKvKIqSS6iCryiKkkuogq8oipJLqIKvKIqSS6iCryiKkkuItKehzRqEEOHAledc3AO4a8E4zysr5MgKGUDleJLK8biskCMrZADzcpSUUnqm9USWLvjmEEKESCkDVY6skUHlUDmyQ46skMGaOdQuHUVRlFxCFXxFUZRcIicX/HlaB0iVFXJkhQygcjxJ5XhcVsiRFTKAlXLk2H34iqIoyuNy8ha+oiiK8ohsX/BFit1CiOaPPNZZCPGXEKKZEOKMEOK8EGK0RhkWCiHuCCFCrdV/JnJsFUJsF0KcEkKcEEIM1SjHX0KIA0KIo6k5JmiVI/V7vRDisBDiD61yCCEuCyGOCyGOCCFCNMqQXwixSghxOvV3pIYGObanvgf/3u4LIYZpkOMvIcTw1N/PUCHET0IIJ41yDE3NcMJi74WUMtvfAF/gFOAE5AHOAS8BF4AXAQfgKFBegwx1gUpAqIbvRS2gUurzrsBZa74XGbwfeVOftwf2A9W1yJH63AjgR+APDX9HLwMeGv5uvAQsBt5MbeMA5NfqZ5L6vB64Rcr55LbOUQu4BDintlkB9NEgRwUgFHAhZRj7LUBpc/vKEePhSylDhRDrgFGkvGFLgELAeSnlRQAhxM9AG+CkrTJIKS8AF4QQ3tbo8xly7Hnk+RghxCmgGFZ6L56S48IjTexTb1Y9iJReDiFEcaAFMImUwm9VT8lh7a6fmgEIJ2WjpE9qmyQgydY5nvjdaAhckFI+70WXz52DlH/AdoCzECKZlIJ7Q4McZYC/pZTxAEKInUA74CtzO8sRt9Q36gxwHHAEOgLfP/J8T+AbW2Z45HFvbLSF/7Qcj2S5CrhpkYOUrbcjQCzwpVbvB7AKqAwEYYMt/KfkuAQcAg4C/WydAQgADgA/AIeB74E8Gv+OLgQGafgzGZr6+xkOLNciB/AKKZ/E3Un5p7MPmGVuPzliCx9AShknhPgFiJVSJoq0N52svTX5WAZr9vU8OYQQeYFfgWFSyvta5JBSGoEAIUR+YLUQwldKadXjG2n8brQE7kgpDwohgqzZ99NypD5cS0p5QwhRCNgshDgtpQy2VQYhhB0puxwHSyn3CyFmAKOBj6yVIa0c/z4uhHAAWgMfWLP/9HIIIQqQsiegFBAFrBRCvCalXGbLHMApIcSXwGZS/vkcBQzm9pPtD9o+wZR6AwgDvB55rjhW/miWRgYtPZZDCGFPSrFfLqX8Tasc/5JSRgE7gGYa5KgFtBZCXAZ+BhoIIaz6B51ODqSUN1K/3gFWA7aY9/DJv5MwKeX+1PurSPkHYAtp/W40Bw5JKW/bKMOTORoBl6SU4VLKZOA3oKYGOZBSLpBSVpJS1gUiSdm3b5acVvAf9Q9QWghRKnWroSuwVuNMmkj9tLMAOCWlnKphDs/ULXuEEM6k/HGdtnUOKeUHUsriUkpvUn4vtkkpX7N1DiFEHiGE67/fA01IOVBnM1LKW8A1IUTZ1IcaYsVjO5nQDfhJw/6vAtWFEC6pfzcNSTmganOpn/oQQpQA2mOB9yXH7NJ5kpTSIIQYBGwkZb/xQinlCVvnEEL8RMp+Yg8hRBgwTkq5wMYxapFyDOO4EOJI6mNjpJTrbZyjCLBYCKEnZWNjhZTSqqdEZnEvkLJbC1L+Fn+UUv6lQY7BwPLUDaOLwOsaZEAI4QI0Bvpr0T9A6m6tVaQcVzGQclxDq6tvfxVCuAPJwEAp5T1zV6iutFUURcklcvIuHUVRFOURquAriqLkEqrgK4qi5BKq4CuKouQSquAriqLkEqrgK4qi5BKq4CuKouQSquAriqLkEv8DTSZ/+ZyTh2kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "outcomes_all_periods = ['Y'+ str(x) for x in periods]\n",
    "\n",
    "plt.plot(outcomes_all_periods, data.loc[treated_unit, outcomes_all_periods].values, label = 'True Untreated')\n",
    "plt.plot(outcomes_all_periods, synthetic_unit[-len(periods):], label = 'Synthetic Control', marker = '.')\n",
    "plt.plot(outcomes_all_periods, data.loc[treated_unit+'treat', outcomes_all_periods].values, label = 'Treated')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "03c99660",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f8de191c1f0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYYklEQVR4nO3df5BV5Z3n8ffHpkkbwDCBNmYEaeJIDIJpsVEZo8mUv0ANaGkKHGcSdmMoJstOkimtQadGmUxS0SFZE0ZYlknYXUczZMsIEiUjFVclUaM02BoREdQY2180bUA6A2Lrd/+4F/Zyvd33dHP719OfV9Ut7znn6ed8Tx/8cHjuOc9VRGBmZgPfUX1dgJmZVYYD3cwsEQ50M7NEONDNzBLhQDczS4QD3cwsEQ50G9QkXS1pfca2cyX9qqdrMusuB7oNSJKul7SuaN32DtbN6aifiLgzIi6sUE0PSbqmEn2ZdYcD3QaqDcDZkqoAJB0HVANTitb9Sb6tWfIc6DZQbSQX4PX55XOBB4FtReteAP4g6UeSXpf0qqRvFYT+YcMoki6UtE3SHknLJD1cfNUt6buSfi/pJUkz8uu+DZwD3CapTdJtyrlV0s58f09LmtRjvxEb9BzoNiBFxAHgcXKhTf6/vwR+VbRuA/C/gXZyV+unARcCHxgakTQauAu4HhhF7i+HPy1qdmZ+/Wjgn4AfSVJE/F1+/wsiYnhELMjv51xgAjASmA20HuGhm3XIgW4D2cP8//A+h1yg/rJo3cPADODrEfGHiNgJ3AqUGle/GNgSEXdHRDuwBHijqM3LEfEvEfEeub8oPg58rIP63gVGACcDioitEfF6N47TLBMHug1kG4DPSPojoDYitgOPAn+aXzcJeI7c0MzrknZL2g38D+DYEv39MfDKwYXIzVzXXNTmjYLt/5F/O7xUcRHxf4HbgKXAm5JWSDqmy0dplpED3Qayx4CPAPOARwAi4m3gtfy614DfAe8AoyNiZP51TEScUqK/14ExBxckqXA5gw9MXRoRSyLidOAUckMv13WhP7MucaDbgBUR+4BG4G/IDbUc9Kv8ug35IY71wPckHSPpKEknSvpsiS7vAyZLukzSEOC/AMd1oaQ3gU8cXJA0VdKZkqqBPwD7gfe60J9ZlzjQbaB7mNzwSeEDP7/Mrzt4u+IXgaHAs8DvyX3w+fHijiJiF/AFch92tgITyf2F8U7GWn4AXJm/A2YJcAzwL/l9vpzv87tdODazLpG/4MKsNElHkRtDvzoiHuzreszK8RW6WQFJF0kaKelDwA2AgF/3cVlmmTjQzQ43jdzDSLuAzwOX5cfqzfo9D7mYmSXCV+hmZokY0lc7Hj16dNTV1fXV7s3MBqRNmzbtiojaUtv6LNDr6upobGzsq92bmQ1Ikl7uaJuHXMzMEpEp0CVNz08pukPSwg7afE5Sk6Qtkh6ubJlmZlZO2SGX/LzRS4ELyD1ksVHS2oh4tqDNSGAZMD0ifiep1MRHZmbWg7KMoZ8B7IiIFwEkrQJmkXuM+qA/B+6OiN8B5Kco7bJ3332X5uZm9u/f350ftw7U1NQwZswYqqur+7oUM+tBWQL9eAqmFCV3lX5mUZsJQLWkh8jN//yDiLi9uCNJ88jNgscJJ5zwgR01NzczYsQI6urqyE10Z0cqImhtbaW5uZnx48f3dTlm1oOyjKGXStbip5GGAKcDlwAXAX8vacIHfihiRUQ0RERDbe0H77rZv38/o0aNcphXkCRGjRrlf/WYDQJZrtCbgbEFy2PIzTNd3GZXRPyB3Pc3bgA+DTzf1YIc5pXn36nZ4JDlCn0jcJKk8ZKGkvvqrrVFbe4BzpE0RNKHyQ3JbK1sqWZm1pmyV+gR0S5pAXA/UAWsjIgtkubnty+PiK2S/h14Gngf+GFEPHOkxdUtvO9IuzjMb2++pNPtra2tnHfeeQC88cYbVFVVcXBo6IknnmDo0KEVrcfMrJIyPSkaEeuAdUXrlhctLwYWV6603jdq1CiampoAWLRoEcOHD+faa689tL29vZ0hQ/rs4VqzfqvSF19HosMLt0Uf6d1COrNoT49063QqY+7cuXz0ox/lySefZMqUKYwYMeKwoJ80aRL33nsvdXV13HHHHSxZsoQDBw5w5plnsmzZMqqqqvr4CMxssPCj/xk8//zz/OIXv+B73/teh222bt3KT37yEx555BGampqoqqrizjvv7MUqzWyw8xV6Bl/4whfKXmk/8MADbNq0ialTpwKwb98+jj3WD8yaWe9xoGcwbNiwQ++HDBnC+++/f2j54P3dEcGXvvQlvvOd7/R6fWZm4CGXLqurq2Pz5s0AbN68mZdeegmA8847j7vuuoudO3OzHrz11lu8/HKHs1yamVVcv75CL3ebYV+44ooruP3226mvr2fq1KlMmJB7IHbixIl861vf4sILL+T999+nurqapUuXMm7cuD6u2MwGi34d6H1p0aJFJdcfffTRrF+/vuS22bNnM3v27B6sysysYx5yMTNLhAPdzCwRDnQzs0Q40M3MEuFANzNLhAPdzCwR/fu2xUrPjpZhhrOqqiomT558aHnOnDksXLiwZNs1a9YwYcIEJk6cCMCNN97Iueeey/nnn39EZe7evZsf//jHfPWrX+3Sz5WaIdLMBo/+Heh94Oijjz40hW45a9as4dJLLz0U6N/85jcrUsPu3btZtmxZlwPdzAY3D7lktHDhQiZOnMipp57Ktddey6OPPsratWu57rrrqK+v54UXXmDu3LncddddQG6KgBtuuIFp06bR0NDA5s2bueiiizjxxBNZvjw3lXxbWxvnnXceU6ZMYfLkydxzzz2H9vXCCy9QX1/PddddB8DixYuZOnUqp556KjfddNOhur797W/zyU9+kvPPP59t27b18m/FzPoTX6EX2bdvH/X19YeWr7/+ei644AJWr17Nc889hyR2797NyJEjmTlzJpdeeilXXnllyb7Gjh3LY489xje+8Q3mzp3LI488wv79+znllFOYP38+NTU1rF69mmOOOYZdu3Zx1llnMXPmTG6++WaeeeaZQ/9SWL9+Pdu3b+eJJ54gIpg5cyYbNmxg2LBhrFq1iieffJL29namTJnC6aef3gu/JTPrjxzoRUoNubS3t1NTU8M111zDJZdcwqWXXpqpr5kzZwIwefJk2traGDFiBCNGjKCmpobdu3czbNgwbrjhBjZs2MBRRx3Fq6++yptvvvmBftavX8/69es57bTTgNyV/fbt29m7dy+XX345H/7whw/bn5kNTh5yyWDIkCE88cQTXHHFFaxZs4bp06dn+rkPfehDABx11FGH3h9cbm9v584776SlpYVNmzbR1NTExz72sUPT8RaKCK6//nqamppoampix44dfPnLXwZAUgWO0MxS4EDPoK2tjT179nDxxRfz/e9//9AV/IgRI9i7d2+3+92zZw/HHnss1dXVPPjgg4em2y3u96KLLmLlypW0tbUB8Oqrr7Jz507OPfdcVq9ezb59+9i7dy8/+9nPun+QZjbg9e8hlx76ItXOFI+hT58+na997WvMmjWL/fv3ExHceuutQO6Wxq985SssWbLk0IehXXH11Vfz+c9/noaGBurr6zn55JOB3JdVn3322UyaNIkZM2awePFitm7dyrRp0wAYPnw4d9xxB1OmTGH27NnU19czbtw4zjnnnCP/BZjZgKWI6JMdNzQ0RGNj42Hrtm7dyqc+9ak+qSd1/t1aT6pbeF9fl3BIh9+jUOnnWo7EEVysStoUEQ2ltnnIxcwsEQ50M7NE9LtA76shoJT5d2o2OPSrQK+pqaG1tdUBVEERQWtrKzU1NX1dipn1sEx3uUiaDvwAqAJ+GBE3F23/HHAP8FJ+1d0R0eWJTcaMGUNzczMtLS1d/VHrRE1NDWPGjOnrMsysh5UNdElVwFLgAqAZ2ChpbUQ8W9T0lxGR7RHKDlRXVzN+/Pgj6cLMbNDKMuRyBrAjIl6MiAPAKmBWz5ZlZmZdlSXQjwdeKVhuzq8rNk3SU5J+LumUUh1JmiepUVKjh1XMzCorS6CXmiyk+FPLzcC4iPg08M/AmlIdRcSKiGiIiIba2touFWpmZp3LEujNwNiC5THAa4UNIuLtiGjLv18HVEsaXbEqzcysrCyBvhE4SdJ4SUOBOcDawgaSjlN+2j9JZ+T7ba10sWZm1rGyd7lERLukBcD95G5bXBkRWyTNz29fDlwJ/JWkdmAfMCd8M7mZWa/KdB96fhhlXdG65QXvbwNuq2xpZmbWFf3qSVEzM+s+B7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmicgU6JKmS9omaYekhZ20myrpPUlXVq5EMzPLomygS6oClgIzgInAVZImdtDuFuD+ShdpZmblZblCPwPYEREvRsQBYBUwq0S7/wr8FNhZwfrMzCyjLIF+PPBKwXJzft0hko4HLgeWd9aRpHmSGiU1trS0dLVWMzPrRJZAV4l1UbT8feBvI+K9zjqKiBUR0RARDbW1tRlLNDOzLIZkaNMMjC1YHgO8VtSmAVglCWA0cLGk9ohYU4kizcysvCyBvhE4SdJ44FVgDvDnhQ0iYvzB95L+F3Cvw9zMrHeVDfSIaJe0gNzdK1XAyojYIml+fnun4+ZmZtY7slyhExHrgHVF60oGeUTMPfKyzMysq/ykqJlZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJcKCbmSXCgW5mlggHuplZIhzoZmaJyBTokqZL2iZph6SFJbbPkvS0pCZJjZI+U/lSzcysM0PKNZBUBSwFLgCagY2S1kbEswXNHgDWRkRIOhX4P8DJPVGwmZmVluUK/QxgR0S8GBEHgFXArMIGEdEWEZFfHAYEZmbWq7IE+vHAKwXLzfl1h5F0uaTngPuA/1yqI0nz8kMyjS0tLd2p18zMOpAl0FVi3QeuwCNidUScDFwG/GOpjiJiRUQ0RERDbW1tlwo1M7POZQn0ZmBswfIY4LWOGkfEBuBESaOPsDYzM+uCLIG+EThJ0nhJQ4E5wNrCBpL+RJLy76cAQ4HWShdrZmYdK3uXS0S0S1oA3A9UASsjYouk+fnty4ErgC9KehfYB8wu+JDUzMx6QdlAB4iIdcC6onXLC97fAtxS2dLMzKwr/KSomVkiHOhmZolwoJuZJcKBbmaWCAe6mVkiHOhmZolwoJuZJcKBbmaWCAe6mVkiHOhmZolwoJuZJcKBbmaWCAe6mVkiHOhmZolwoJuZJcKBbmaWCAe6mVkiHOhmZolwoJuZJSLTd4r2N3UL7+vrEg757c2XdLxx0Ud6r5DOLNrT1xWYWS/wFbqZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmicgU6JKmS9omaYekhSW2Xy3p6fzrUUmfrnypZmbWmbKBLqkKWArMACYCV0maWNTsJeCzEXEq8I/AikoXamZmnctyhX4GsCMiXoyIA8AqYFZhg4h4NCJ+n1/8NTCmsmWamVk5WQL9eOCVguXm/LqOfBn4eakNkuZJapTU2NLSkr1KMzMrK0ugq8S6KNlQ+jNygf63pbZHxIqIaIiIhtra2uxVmplZWVnmcmkGxhYsjwFeK24k6VTgh8CMiGitTHlmZpZVliv0jcBJksZLGgrMAdYWNpB0AnA38JcR8XzlyzQzs3LKXqFHRLukBcD9QBWwMiK2SJqf374cuBEYBSyTBNAeEQ09V7aZmRXLNH1uRKwD1hWtW17w/hrgmsqWZmZmXTEg50O3dPSXue09r72lwI/+m5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klIlOgS5ouaZukHZIWlth+sqTHJL0j6drKl2lmZuUMKddAUhWwFLgAaAY2SlobEc8WNHsL+Gvgsp4o0szMystyhX4GsCMiXoyIA8AqYFZhg4jYGREbgXd7oEYzM8sgS6AfD7xSsNycX2dmZv1IlkBXiXXRnZ1JmiepUVJjS0tLd7owM7MOZAn0ZmBswfIY4LXu7CwiVkREQ0Q01NbWdqcLMzPrQJZA3wicJGm8pKHAHGBtz5ZlZmZdVfYul4hol7QAuB+oAlZGxBZJ8/Pbl0s6DmgEjgHel/R1YGJEvN1zpZuZWaGygQ4QEeuAdUXrlhe8f4PcUIyZmfURPylqZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmiXCgm5klwoFuZpYIB7qZWSIc6GZmicgU6JKmS9omaYekhSW2S9KS/PanJU2pfKlmZtaZsoEuqQpYCswAJgJXSZpY1GwGcFL+NQ/47xWu08zMyshyhX4GsCMiXoyIA8AqYFZRm1nA7ZHza2CkpI9XuFYzM+vEkAxtjgdeKVhuBs7M0OZ44PXCRpLmkbuCB2iTtK1L1fZDuuWwxdHArr6ppBP/oEr21j+P8QgVnUfoj8dZ2fMI/fEYK2AQ/D85rqMNWQK91J6jG22IiBXAigz7HJAkNUZEQ1/X0ZMGwzHC4DhOH2N6sgy5NANjC5bHAK91o42ZmfWgLIG+EThJ0nhJQ4E5wNqiNmuBL+bvdjkL2BMRrxd3ZGZmPafskEtEtEtaANwPVAErI2KLpPn57cuBdcDFwA7gP4D/1HMl92vJDicVGAzHCIPjOH2MiVHEB4a6zcxsAPKTomZmiXCgm5klwoHeDYNhKoQMx/g5SXskNeVfN/ZFnUdC0kpJOyU908H2FM5juWNM4TyOlfSgpK2Stkj6Wok2A/5cZhIRfnXhRe6D4ReATwBDgaeAiUVtLgZ+Tu7+/LOAx/u67h44xs8B9/Z1rUd4nOcCU4BnOtg+oM9jxmNM4Tx+HJiSfz8CeD61/yezvnyF3nWDYSqELMc44EXEBuCtTpoM9POY5RgHvIh4PSI259/vBbaSe1K90IA/l1k40Luuo2kOutqmP8ta/zRJT0n6uaRTeqe0XjXQz2NWyZxHSXXAacDjRZsGxbnM8ui/Ha5iUyH0Y1nq3wyMi4g2SRcDa8jNtpmSgX4es0jmPEoaDvwU+HpEvF28ucSPpHYufYXeDYNhKoSy9UfE2xHRln+/DqiWNLr3SuwVA/08lpXKeZRUTS7M74yIu0s0Sf5cggO9OwbDVAhlj1HScZKUf38GuT9Lrb1eac8a6OexrBTOY77+HwFbI+K/ddAs+XMJHnLpshgEUyFkPMYrgb+S1A7sA+ZE/naCgULSv5G7y2O0pGbgJqAa0jiPkOkYB/x5BM4G/hL4jaSm/LobgBMgnXOZhR/9NzNLhIdczMwS4UA3M0uEA93MLBEOdDOzRDjQzcwS4UC3AUXSe/lZAbfkH1f/G0m9+udY0kOSGgqW6zqazbDo59ZJGpl/fbVnq7TByIFuA82+iKiPiFOAC8jdW3xTT+5QUlUl+omIiyNiNzAScKBbxTnQbcCKiJ3APGBB/gnAGkn/U9JvJD0p6c8AJM2VdLekf5e0XdI/HexD0lX59s9IuqVgfZukb0p6HJiWtaYy+/pt/rH6m4ET8//SWFyBX4UZ4CdFbYCLiBfzQy7HAn+RXzdZ0snAekkT8k3ryc3C9w6wTdI/A+8BtwCnA7/Pt78sItYAw8jNId6dL3z4wL4ionCmv4XApIio70bfZh3yFbql4OBMep8B/hUgIp4DXgYOBvoDEbEnIvYDzwLjgKnAQxHREhHtwJ3kvhACcmH/0w72V+rx6sJ1pfZl1uMc6DagSfoEufDdSekpUg96p+D9e+T+ddpZ+/0R8V4H21qBPypY/iiwq8y+zHqcA90GLEm1wHLgtvyEUhuAq/PbJpCbnGlbJ108DnxW0uj8B59XAQ9n2PVDwF8cnKUQ+BLwYBdK30vuq9LMKsqBbgPN0QdvWwR+AawH/iG/bRlQJek3wE+AuRHxTgf9kJ8+9XpyYfwUsDki7slQwwpyofyUpKeA4cB3sx5ARLQCj+Q/iPWHolYxnm3RzCwRvkI3M0uEA93MLBEOdDOzRDjQzcwS4UA3M0uEA93MLBEOdDOzRPw/pbpSG7jJGf0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(np.arange(len(true_weights))-0.18, true_weights, width = 0.3, label = \"True\")\n",
    "plt.bar(np.arange(len(true_weights))+0.18, estimated_weights, width = 0.3, color = 'C1', label = 'Estimated')\n",
    "plt.title(\"Weights\")\n",
    "plt.xlabel(\"Donor Unit\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7043b77f",
   "metadata": {},
   "source": [
    "## Calculate $p$-values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2ea676f8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/_r/5whyybsj5jlg0yywjfk54qpr0000gp/T/ipykernel_30489/471669146.py:2: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "  ratios = pd.Series()\n"
     ]
    }
   ],
   "source": [
    "placebo_effects = dict()\n",
    "ratios = pd.Series()\n",
    "\n",
    "treated_unit_treat = treated_unit + 'treat'\n",
    "for unit in list(donor_units) + [treated_unit_treat]:\n",
    "\n",
    "    synth_outcomes, true_outcomes, V_star, estimated_weights = synthetic_control(data = data,\n",
    "                            features = features,\n",
    "                            treated_unit = unit,\n",
    "                            units = [x for x in list(donor_units)+[treated_unit] if x!= unit], # donors\n",
    "                            periods = periods,\n",
    "                            pretreatment_periods = pretreatment_periods,\n",
    "                            n_features = n_features,\n",
    "                            search_size = 100)\n",
    "        \n",
    "    placebo_effects[unit] = synth_outcomes, true_outcomes\n",
    "    ratio = rmspe_ratio(synth_outcomes, true_outcomes, treatment_start = 7)\n",
    "    ratios.loc[unit] = ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "b5ef11f8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.25 1 4\n"
     ]
    }
   ],
   "source": [
    "ratios.sort_values(ascending = False, inplace = True)\n",
    "\n",
    "rank = list(ratios.index).index(treated_unit_treat) + 1 \n",
    "\n",
    "p_value = rank / len(ratios)\n",
    "print(p_value, rank, len(ratios))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06dde06f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}