1.5 - Structural Causal Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"width: 100%; overflow: hidden;\">\n",
    "    <div style=\"width: 150px; float: left;\"> <img src=\"https://raw.githubusercontent.com/DataForScience/Networks/master/data/D4Sci_logo_ball.png\" alt=\"Data For Science, Inc\" align=\"left\" border=\"0\" width=150px> </div>\n",
    "    <div style=\"float: left; margin-left: 10px;\"> <h1>Causal Inference In Statistics - A Primer</h1>\n",
    "        <h1>1.5 Structural Causal Models</h1>\n",
    "        <p>Bruno Gonçalves<br/>\n",
    "        <a href=\"http://www.data4sci.com/\">www.data4sci.com</a><br/>\n",
    "            @bgoncalves, @data4sci</p></div>\n",
    "    <div style=\"float: right; margin-right:10px;\"> <p><a href=\"https://amzn.to/3gsFlkO\" target=_blank><img src='data/causality.jpeg' width='100px'>\n",
    "        <!--Amazon Affiliate Link--></a></p></div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "from pprint import pprint\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "import statsmodels\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "from CausalModel import CausalModel\n",
    "\n",
    "import watermark\n",
    "\n",
    "%load_ext watermark\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by print out the versions of the libraries we're using for future reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy           1.18.1\n",
      "watermark       2.0.2\n",
      "statsmodels     0.11.1\n",
      "matplotlib      3.1.3\n",
      "autopep8        1.5\n",
      "json            2.0.9\n",
      "statsmodels.api 0.11.1\n",
      "pandas          1.0.1\n",
      "Sat Jul 25 2020 \n",
      "\n",
      "CPython 3.7.3\n",
      "IPython 6.2.1\n",
      "\n",
      "compiler   : Clang 4.0.1 (tags/RELEASE_401/final)\n",
      "system     : Darwin\n",
      "release    : 19.5.0\n",
      "machine    : x86_64\n",
      "processor  : i386\n",
      "CPU cores  : 8\n",
      "interpreter: 64bit\n",
      "Git hash   : 8e58bf1ca04167a82c7ab09e68e7b88b5ae50adb\n"
     ]
    }
   ],
   "source": [
    "%watermark -n -v -m -g -iv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load default figure style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('./d4sci.mplstyle')\n",
    "colors = plt.rcParams['axes.prop_cycle'].by_key()['color']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modeling Causal Assumptions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can load the structure of this model directly from disk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = CausalModel('dags/Primer.Fig.1.9.dot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot it quickly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAugAAAH4CAYAAAAYZBiWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAuIwAALiMBeKU/dgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeVzUdeIG8GeGQ0RURC1NDe/7YjgUD7zPzLQZDu8DzStNK680S3NXLcstWy1bzVqPfu7WttvaWnlFCIIwaOaNeKMoct8w8/394TbbJ1ARYT4z33ner9e+9tUjjg+H8jDzmflqFEVRQERERERENkEruwAREREREf0PBzoRERERkQ3hQCciIiIisiEc6ERERERENoQDnYiIiIjIhnCgExERERHZEA50IiIiIiIbwoFORERERGRDONCJiIiIiGwIBzoRERERkQ3hQCciIiIisiEc6ERERERENoQDnYiIiIjIhnCgExERERHZEA50IiIiIiIbwoFORERERGRDONCJiIiIiGwIBzoRERERkQ3hQCciIiIisiEc6ERERERENoQDnYiIiIjIhnCgExERERHZEA50IiIiIiIbwoFORERERGRDONCJiIiIiGwIBzoRERERkQ3hQCciIiIisiEc6ERERERENoQDnYiIiIjIhnCgExERERHZEGfZBdQkr1jBpfR7/8spUlBkAopMgEkBnLVANSfAxUmDutWBFl4aNKmtgbNWI7s2ERERkaDErOBapoKLaQru5gPFJgWFJqDEDDhpAFene//zcNWgWZ17/3N34aapLBzoFZRfrCAp/d4X7sV0BUlpZiRnA8oj3IaLFvD21KCFlwYt6tz7/6dra+DE0U5ERERWYjIruJr5v01zMU3BlQwFxeby34YGwFM1geZeWsumaV5Hg+oc7RWiURTlUTalw7tw14x9F8yIvGp+pC/c8qrjBgxqocWgFk6o684vaiIiIqoaqXkK9l804YeLZqQXVP7tuzgBvZ7WYmhLLVrV5anqR8GBXg6FJQoir5qxL9GMi2nW+XBpNUBAIw2GtXJCxyc00Gg41omIiOjxKIqCkykK9iWaEHtDgdlKK7CFlwZDW2rR62ktqjlz0zwMB/oD3M1T8K9zJhy6ZEZOkbwejWoCQ1s5YXALLVyc+EVNREREj6bYpOD7i2bsu2DCjWx5PTxcgX7NtBjZhicFHoQDvQyKouBAkhmfHjchv1h2m/9pUkuDF7s58WEiIiIiKrcLd834MMaEa1m2M/mquwBTfZzQv5mWpwTKwIH+O6m5CjbHlSDhpm1+WLQaYFRbLUI7OvHedCIiIrqvYpOC//vFhH+cMT/Si1hYk09DDWb5O6Me700XcKD/l63ea34/vDediIiI7scW7zW/H96bXhoHOoC0fAV/jrXde83v59d708d0cuJLMxIRERFMZgW7T9r2veb349NQgzkBzvCqzk3j8AP9VraClYeLkZIru0nFBTTS4OUeznDlkRciIiKHVWRS8F5UCWJv2O+0e7IG8EY/FzTwcOxN49AD/UqGGSsPlyCjCl7709o6P6nB4l7OvCAAERGRA8ovVrD2pxKcvG3/s66OG7CirzO8PR33GK/DDvTLGWasOFgi9eUTK1ubuhqs6MuRTkRE5EjyixWsOlyCc3fVM+k8XIFV/Z3R1EFHukO+18nZ976Q1TTOAeDcXQVvR5ag2KSev6BERER0f8UmBesi1TXOASCnCFh1uAQ3s9X1fpWXww301DwFKw8Vq+JYS1lOpCh4L7oEJmtdGoyIiIikMJnvfc//OUWd3/MzCoA3DxUjNU+d79+DONRAN5kVvBtVgjt5sptUrZjrCr46Y5Zdg4iIiKrQl6fNiLmu7vF6Jw94L8rx7nh0qIG+97wZ51Id4xO85xcTLmdwpBMREanR5Qwz/nbKJLuGVZxNVfDtBcfaNA7zJNHkLAUL9hWj+DE/v2f+tgpnv/qDkNVp1Q193jwEjdapwrerKAqi1jyD2ycPCHmLoXPQedJ7FbrN5nU0WDvIGc58jXQiIiLVKDErWPxDCS6lV2zCKWYTflzRB+kXjwl5W/3raGdYXhkVUZh1BwcW6lCYdVvIuy/8Cg11zzzy7blogQ3DXPBUTcfYNA5xD7rJrGBjbMljj3MAaDN6KWp7dxGy9AsxuPDvDY91u5d++LjUOPdo2AodxvzhPr/j4ZLSFfyDR12IiIhU5avT5gqPcwDQaJ3gO3srnFyrC/m5r9ci41LC49YDABzfOrfUOPfuO7lC4xwAis3AhzGOc9TFIQb6txcq72iL1tkFvrO2QuvsKuRn/r4KWddOV+g2c1Iu4pddrwmZRusE31ml//I8qj2/mHCFR12IiIhUobKOttR8qg3ah70lZIqpGHGbpsJUXPhYt30tcjeSY/8hZO71vNFpwjuPdbuOdNRF9QM9OUvBjhOVe0artncntNWLDwGZiwsRvzkcZlPJI92WYjYjfvM0mArFS5m2HLEAXq26PXZXkwJsjDGhxEF+4iQiIlKrErOCD2NMqKxXU24x9EXUa99HyLKvn8bpPW9U+Dbz05JxYvsCMdRooJu5BS7utSp8u7/a+bMJyQ7w0ouqH+i7fzFVytGW32s98lXUaeEvZBmXjDj/9bpHup0Lezcg7VyUkNVq0gHtgyv+l+P3ktIVRF5xjJ84iYiI1OqnK2YkPcbRlt/TaDTwnfkJnKvXFPLEve8j9eyRCt1mwiczUZybLmQthsxB/Q59K1pTUGQCvjip/ifHqnqgp+crOHqtaoapRusE3axPoHVxE/Kz/1iDjMsnynUbWddO48zfVoq361T2EZrHtS+RA52IiMieVcX3cvf63ug0Yb0YKmYYN4ejpCDnkW7r0oGtSDn+nZB5PNUaHcasftyaguhrZqTnq/tedFUP9P1J5kp7GKgstRq1Q/uQN4VMMRUjftNUmEsefJlSs6kE8Zunwvy7c15tRi2BZzOfyq6K83cVXEzjSCciIrJHF9PMuFBFVwtt2m8yGuiGC1nu7Us4uWNxuW8j985l/LJTfPt7z6fb9tjPp/s9kwIcSFL3plHtQDeZFXyfWPUPgbQc/hK82vQQsqxrv+DM39+6z++459w/1pR6prRnMx+0Gb2k0jv+6jvei05ERGSXqvqRcJ/pm+HqUVfILh/4S6l7xMuiKAqMH72AkvxsIW/93CJ4tfS/z+96PN9fNKn6FV1UO9DjkhXcza/6P0ej1cJ35idwquYu5Be+eRdpicfK/D3pSUac+3qtkGldqsF31jZonZyrrGvEFTNyi9T7xUxERKRGOUUKfqri55K5eTZA1/CNpXLjlpkoykkv43f8z8V9HyL19I9CVtu7C9o+v6xSO/5Wah4Qn6zeTaPagb7vgvWeQODRoGWp1ytXzCbEbw6HqahAyE3FhYjfNBXK717tpZ1hBWo1aV+lPYtMwKFLvBediIjInhy6ZEaRFWZNo+56NO4RKmQF6ck4sX3+fX9P9s3zOP3F60KmdakGvzmfQuvsUiU9f7XPCiclZFHlQL+VreBEinV/qmo+eBbqte8rZDnJ53D6/1YI2ek9byD7xhkh82rVHa1G/O4liarIdxfV+8VMRESkRtY4svurLlPeh1udp4Ts+pEvcCPmq1Jve+/OyGkwFYlHFu7d6dihSnsCwPFbCm6p9CUXVTnQf7lt/XuJ771U0ZbSL1X0n42WlypKPXsEiXvfF37dqZo7fGf9BRqtk1V63siC6p/5TEREpBZp+QpuZD/87SqLq0cd+LywuVR+fOtcFGSkCNn5b95D+oUYIfNq0wOtRrxcpR1/69QddZ4MUOVAv1iJrxH6KNzre6PjuN+9DrpihvGjaSjMvA3jR9MARfxC6hC2Gh4NW1mxpbyPDxERET2apDTrf89u0HUomg4IF7Ki7FQkfDLL8t9Z107h7N9XCW/jVK3Gf+90tN68vCjh42MNqhzoMr6Yf9VsQDie6DxIyHJTknBgiT9yU5KEvF77vmg+ZLY16wGQ+/EhIiKi8pN1p1qn8W/D/YmmQnbLuBdXDn8Gc0kx4sp4SemO49bA48kWVmyp3jsdVTfQS8wKLmfI/WTpZnwMlxqeQlaYcUv4b+fqNeE7cws0Go01qwEAEvl66ERERHZB1jVMnN084DvzL4BGnIo/f/4qTnw6D5mXjwv5E50HovmgGdasCAC4nK6gRIUvt6i6gX4tU0Gx5P1Z3asROk9894Fv03HcOrjX97ZSI1FlXiaYiIiIqo7MIxz12vVGy+HzhKwkPwuXD24TMpcantC98LE1q1kUm+9tP7VR3UC3lbNITweNxxOdB5b5a092HYJmvzvbZU1p+XyiKBERka1Ly1eQXvDwt6tK7UNWomajdg98m86TNqB63cZWalSaGu94VN1At5VPUmFWKjKvnCzz154KeN7KbUpT65ktIiIitbCFTePk6gbf2duguc+FFJ8KGIWne4+1ciuRrdw5W5lUN9AzCmzjk3Ri2zwUZqaU+WunvliOwqw7Vm4kyrSRjxMRERGVzVa+V9dprkObUUtK5a4166Fr+IcSGokyC23j41SZVDfQrXGlrYe5HrUHN2K+vO+vF2XdwYltL1mxURkdbODjRERERPdnS9+r67UPKpXVatwe1WrVl9BGZEsfp8qiuoFeLPmTVJBxC8c/LT2+XdxrC/99I+ZL3Dh6/xFf1WR/nIiIiOjB1Dg8q4IaN436BrrkV3BJ+GQ2inPShKzZwOnQzfyk1Nse//QlaUddZH+ciIiI6MHUODyrgho/Tqob6M4S36MrP36OW8a9Qub+RFN0HLcWT/k/h8Y9QoVfK8q6g+PbxJcvshaZHyciIiJ6OGcn2Q3sgxo3jereJVdJX8z5d6/j5OeviqFGA98Zn8DZzQMA0GXKn1DNs4HwJskxX+F69N+tVdNC1seJiIiIyqcav1eXi2vZLzBj1zjQK4lxywwU52UKWYshs4UnVbh6eKFr+MZSv/fEpy+hMPN2lXf8LQ50IiIi2+aiupVWNdT4cVLdu1Srmsbqf+alH7bg9s/7haxGg5ZoH7a61Ns+5TcSTXqJrxdalJ1q9aMuMj5OREREVH613Pi9ujxqq/DjpLqB3ryOdT9Jubcv4ZddS8VQo4XvrL/AuZp7mb+n86T34ObZUMiSY/+B69F/q6qapVj740RERESPht+ry0eNHyfVDfQWXtb7JCmKAuNHL6CkIEfIWz7zEuq2Drzv73P1qIOu00q/sP+JT+db5aiLpxvgVb3K/xgiIiJ6DHWrA7WryW5h+zjQ7cDTtTVWezbvxf9sROqZCCGr2agt2ge/+dDf29B3BJr0Hi9k9466zK3MimVqXkcDjUZ9X8xERERqotForHrHoz1y1t7bfmqjuoHu4qSBtxU+Udk3z+P0/60QMo3WCb6ztsLJ1a1ct9F50rtwq/OUkCXHfo3rUXsqrWdZWvIvOxERkV3gQH8wb08NXJzU9zFS3UAHgOZV/MWsmE2I3zwNpqJ8IW818lXUaeFX7ttxreEJn2l/LpWf+HQ+CjJSHrvn/TT3UuWnnYiISHVa1OH37AdpocLjLYBKB3pVf7LOf/Me0i/ECFmtJh3RTr/8kW+rgW44nu4zUciKcu5W6VEXtX4xExERqU1V3+lo79T6CIMqB3qHJ6ru3cq6dgpn/75KyDROLvCdtRVaZ9cK3WanCe/AzauRkN089k9cO/J/Fe55Pw08+ARRIiIie1G3OvBkDdktbFf7+qqcsuoc6I1qadChfuX/RGUuKUb85nCYS4qEvM2oJfBs1rXCt+tawxM+0zeVyn/evqDSj7oMbuHEJ4gSERHZCY1Gg8EteXXBsnR8QoNGtdS5aVQ50AFgaKvKf9fOfb0WGZcShMyzmQ/ajF7y2LfdoOtQePedLGRFOXdxfOuLj33bv3LRAv2bq/ZTTkREpEoDmmtVebXMxzVUxT+4qPbTHdBIC8/yvZhKuWRcSsC5r9cKmdbZ9d7RFifnSvkzOk14G9W9GgvZzbh/4dqRLyrl9ns+reUVRImIiOxMrWoa9HhatZOtQjzdgIDG6t00GkVRFNklqsrukyX42ymz7Bo2Y+0gZ7Suy7/gRERE9uZcqhlL95fIrmEzQjpoEdapcu4gtUWqXmuDWjhBq94frh5J8zoatFLpM52JiIjUrnVdDZrxVdgAAFoNMLCFeo+3ACof6PXcNfBvxC9mABjSUssnhxIREdkpjUaDoS1VPdvKzb+RBvXc1b1pVP+ZDuvoBBVeYOqRNKmtQd+mqv9UExERqVrfplo0UemrlpSXs/betlM71a82b08tQh3gE3k/Wg0wt5uTKi+DS0RE5EhcnDR4sZtjH98N6eAEb0/Vz1f1D3QAGN1O67BXzxzVVouWXg7xaSYiIlK9VnW1GNXWMb+vt6ijweh2jvG+O8R76aS99xOno92J3KS2xqEfPSAiIlKj0I5ODnfUxVmLe1vOQR4+cIiBDjjeURcebSEiIlInRzzq4ihHW37lOO8pHOuoC4+2EBERqZcjHXVxpKMtv3Ko99ZJq8ErPZwr9QqjtqhrAx5tISIiUrvQjk7o2kDddzzWcQNe6ensMEdbfuVQAx0AGtTU4I2+zvBwld2karStp8GiXs482kJERKRyLk73vue3qafO7/kersCKvs5o4KHO9+9BHG6gA/fOoy/v4ww3lV0htpmnBsuCnOHm7HhfyERERI7IzVmD5UHOaOapru/9bs7A8j7ODnXu/Lc0iqIoskvIcuGuGW/9WIKcItlNKk5RFGg0GrSpq8GyPs7wcFXXX1AiIiJ6uOxCBX+IKMH5u4plG9grD1fg9T7OaFXXMcc54KD3oP+qVV0tVg9wRh07PpOu0Wjgnp2IN/pxnBMRETmqmtXuHeF1z75o1+PcqzqweoBjj3PAwQc6ADxdW4s1A13Q1k7Pb138bjN2z/bB999+I7sKERERSfT9t//C7tldcfG7zbKrVEjbehr8cYALnq7t8PPUsY+4/JbJrODbC2bsOGFCsVl2m4fLvX0Jxo9nIPX0jwCA2rVrIz4+Hi1atJDcjIiIiKzt4sWL8PX1RWZmJgCgXvs+8J25Be71m8otVg6uTsC4zk4Y3krrcK/Wcj/8EeW/nLQaPNvGCRuG2f696cNaaTFMe8gyzgEgMzMTBoMB+fn5EpsRERGRteXn50Ov11vGOQCknv4Rw7SHMLSlbU+9tvU0eG+oC55t4zhXCS0P3oNeBlu9N/3JGsDsAGd0elILRVEwduxYfPHFF8LbhIeH4y9/+YukhkRERGRt4eHh2LZtm5CNGTMGO3fuhEajwckUM/4cW4LbuZIKloH3mj8YB/oDJGcr+PqMCRFXzCgyyetRtzowuKUTRrTWorrL/76Ic3JyEBAQgDNnzghvv3XrVkydOtXaNYmIiMjKtm7dimnTpglZu3btEBsbCw8PD0uWX6zgm3NmfH/RhDSJD7a7OgFB3lqMbueEhjU5zO+HA70ccooUHL5kxr5EE5KzrffndnlSg6GtnOD3lOa+P12eOXMG/v7+yM3934/Fbm5uiI6ORteuXa1VlYiIiKwsISEBgYGBKCwstGQ1atTAsWPH0K5duzJ/j8ms4Fiygn0XTPg5xXoT8KmawNCWTujXTIsafNW5h+JAfwSKouBkioL/JJoQe11BVXzgargA/ZtrMaSFE56qVb4v4C+++AJjxowRsubNmyM+Ph6enp5V0JKIiIhkysjIgK+vL5KSkoT8iy++QGhoaLlu40aWgu8STTh4yYy84srvqNUAAY00GNrSCZ2e1Nj1yz9aGwd6BWUUKLhwV8HFNDOS0hUkpinIKHi029Dg3k+Uzb20aOmlQfM6GrT00qBaBa4EOm/ePGzcuFHInnvuOfzjH//gXwgiIiIVMZvNGD16NP71r38J+bx58/D+++8/8u0VltzbMb/umaQ0M5Kz8ch3RHq6wbJnWnhp0aquBp5u3CAVwYFeidLyFVxMU3ApXUFOkYIiE1BsAkrMgIvTvXNXLk5APfd7X7zN62iEM+WPo6ioCH369MHRo0eFfN26dVi0aFGl/BlEREQk37p167BkyRIhCwwMxOHDh+Hq6lopf0Z+8b3BnpSuIDVPQbEJll3jrP3frvFw1aBZHQ1aeGngVZ1jvLJwoKvItWvXoNPpkJqaasm0Wi0OHDiAvn37yitGREREleLw4cMYMGAAzOb/vcxcvXr1kJCQgMaNG0tsRpXJtl8ckx5JkyZNsGvXLuFIi9lsRlhYGG7evCmxGRERET2u5ORkhIaGCuNco9Fg165dHOcqw4GuMoMGDcKbb74pZCkpKQgNDUVxcRU8A4SIiIiqXHFxMUJDQ3H79m0hX7lyJQYNGiSpFVUVDnQVWr58OYYNGyZkP/30E1577TVJjYiIiOhxLF26FJGRkUI2bNgwLFu2TFIjqko8g65Sd+/ehU6nw9WrV4X8q6++wujRoyW1IiIiokf11VdfQa/XC5m3tzeMRiO8vLwktaKqxIGuYseOHUOvXr1QVFRkyWrVqoW4uDi0atVKYjMiIiIqjwsXLsDPzw9ZWVmWzNXVFUeOHIGfn5/EZlSVeMRFxfz9/fGnP/1JyLKysmAwGJCXlyepFREREZVHXl4e9Hq9MM4B4P333+c4VzkOdJWbOXMmxo0bJ2Q///wzZs+eDT54QkREZJsURcGsWbNw8uRJIR8/fjxmzJghqRVZC4+4OIDc3Fx069YNp06dEvItW7Zg+vTpkloRERHR/WzZsqXUEO/YsSOOHj2KGjVqSGpF1sKB7iDOnTsHf39/ZGdnW7Jq1aohKioKOp1OYjMiIiL6rfj4ePTo0UN4DlnNmjURFxeH1q1bS2xG1sIjLg6iTZs22LZtm5AVFhZCr9cjPT1dUisiIiL6rbS0NBgMBmGcA8C2bds4zh0IB7oDMRgMmD9/vpBdvnwZEydOFK5KRkRERNZnNpsxceJEXL58WcgXLFgAg8EgpxRJwYHuYN5++2306NFDyP79739j7dq1khoRERERAKxZswZ79+4Vsp49e2LdunWSGpEsPIPugG7cuAEfHx/cuXPHkmm1Wnz//fcYMGCAxGZERESO6cCBAxg8eLDwiPYTTzwBo9GIRo0aSWxGMnCgO6iDBw9i0KBBwj8E9evXR0JCAv8hICIisqL73XG2f/9+9OvXT2IzkoVHXBxU//798dZbbwnZnTt3EBoaiuLiYkmtiIiIHEtxcTFCQkKEcQ4Aq1ev5jh3YBzoDmzJkiUYMWKEkB05cgSLFy+W1IiIiMixLFq0CFFRUUI2YsQIfi92cDzi4uDS09Oh0+lKPWN8z549CA4OllOKiIjIAezZswehoaFC1qxZM8THx6NOnTqSWpEt4EAnxMfHo2fPnigsLLRkHh4eiIuLQ5s2bSQ2IyIiUqezZ8/C398fOTk5lowXEKRf8YgLwdfXFxs3bhSynJwc6PV65ObmSmpFRESkTrm5uTAYDMI4B4APP/yQ45wAcKDTf02bNg2TJk0SslOnTmHmzJnggyxERESVQ1EUzJgxA6dOnRLyyZMnIzw8XFIrsjU84kIWeXl5CAwMxM8//yzkmzdvxsyZMyW1IiIiUo/Nmzdj9uzZQta5c2dER0fD3d1dUiuyNRzoJLhw4QL8/PyQlZVlyVxdXREZGQl/f3+JzYiIiOxbbGwsevXqJbycca1atRAfH4+WLVtKbEa2hkdcSNCqVSts375dyIqKimAwGHD37l05pYiIiOzc3bt3ERwcXOpaI5999hnHOZXCgU6ljB49Gq+++qqQXb16FRMmTBCuPEpEREQPZzabMX78eFy9elXIFy5ciFGjRklqRbaMA53KtGbNGvTu3VvI/vOf/+APf/iDpEZERET2afXq1di3b5+QBQUF4Y9//KOkRmTreAad7uvmzZvw8fFBSkqKJdNoNPjuu+8waNAgic2IiIjsw/fff4+hQ4cKr4j25JNPIiEhAQ0bNpTYjGwZBzo90OHDhzFgwADhaEvdunWRkJCAJk2aSGxGRERk265evQqdTic8h8vJyQkHDhxAnz59JDYjW8cjLvRAffv2xZo1a4Ts1ye6FBUVSWpFRERk24qKihASElLqBRbWrFnDcU4PxYFOD7Vw4UI899xzQhYTE1PqiaRERER0zyuvvIKYmBghGzVqFL93UrnwiAuVS0ZGBvz8/HDx4kUh3717N8LCwiS1IiIisj1ffPEFxowZI2QtWrRAXFwcPD09JbUie8KBTuV2/PhxBAYGoqCgwJLVqFEDx44dQ7t27SQ2IyIisg2nT59GQEAAcnNzLZmbmxuOHj2KLl26SGxG9oRHXKjcunbtij//+c9ClpubC71ej5ycHEmtiIiIbEN2djb0er0wzgFg06ZNHOf0SDjQ6ZFMnToV4eHhQnbmzBlMnz4dfDCGiIgclaIomD59Os6ePSvk06ZNw5QpUyS1InvFIy70yPLz89GjRw8cP35cyDdu3IgXX3xRUisiIiJ5Nm7ciHnz5gmZj48PoqKi4ObmJqkV2SsOdKqQpKQk6HQ6ZGZmWjIXFxdERESge/fuEpsRERFZ19GjRxEUFITi4mJL5unpifj4eDRv3lxiM7JXPOJCFdK8eXN8/vnnQlZcXIzg4GDcuXNHUisiIiLrunPnDoKDg4VxDgCff/45xzlVGAc6VdjIkSOxePFiIbt+/TrGjRsHk8kkqRUREZF1mEwmjB07FtevXxfyJUuW4Nlnn5XUitSAA50ey+rVq9G3b18h++GHH7Bq1So5hYiIiKxk5cqV2L9/v5D169cPb731lqRGpBY8g06P7datW9DpdLh586Yl02g0+PbbbzF06FCJzYiIiKrGf/7zHwwfPlzIGjZsiISEBDz55JOSWpFacKBTpfjpp5/Qr18/4WiLl5cXjEYjvL29JTYjIiKqXFeuXIFOp0NaWpolc3JywuHDh9GrVy+JzUgteMSFKkXv3r2xbt06IUtLS0NwcDAKCwsltSIiIqpchYWFMBgMwjgHgLfffpvjnCoNBzpVmpdffhnPP/+8kB07dgwLFiyQ1IiIiKhyzZ8/H3FxcUKm1+v5vY4qFY+4UKXKzMyEn58fEhMThXzHjh0YN26cpFZERESPb8eOHZgwYYKQtWrVCnFxcahVq5akVqRGHKTM5kUAACAASURBVOhU6X7++Wd0794d+fn5lszd3R2xsbHo0KGDxGZEREQV88svvyAgIED43la9enXExMSgU6dOEpuRGvGIC1W6zp0746OPPhKyvLw86PV6ZGdnS2pFRERUMVlZWdDr9cI4B4CPP/6Y45yqBAc6VYmJEyfihRdeELJz584hPDwcfNCGiIjshaIoCA8Px/nz54V8xowZpY67EFUWHnGhKlNQUICePXvCaDQK+YYNGzB//nxJrYiIiMpvw4YNePnll4XM19cXkZGRcHNzk9SK1I4DnarUpUuX4Ovri/T0dEvm7OyMw4cPo2fPnhKbERERPdiRI0fQt29flJSUWLI6derAaDSiadOm8oqR6vGIC1WpZs2a4a9//auQlZSUICQkBLdv35bUioiI6MFu376NkJAQYZwD917JheOcqhoHOlW5Z555BsuWLROy5ORkjB07VrjyKBERkS0wmUwYM2YMkpOThXz58uUYPny4pFbkSDjQySpWrlyJAQMGCNmBAwfwxhtvSGpERERUthUrVuDgwYNCNnDgQLz55ptyCpHD4Rl0sprbt29Dp9Phxo0bQv7vf/8bzzzzjKRWRERE//Pvf/8bzz77rJA1atQICQkJqF+/vqRW5Gg40MmqoqKi0KdPH+FMn6enJ4xGI5o1ayaxGRERObqkpCT4+voiIyPDkjk7OyMiIgKBgYESm5Gj4REXsqoePXpg/fr1QpaRkQGDwYCCggJJrYiIyNEVFBQgODhYGOcA8O6773Kck9VxoJPVzZs3DyEhIUJmNBrx0ksvSWpERESObt68eaWu2xESEoK5c+dKakSOjEdcSIrs7Gz4+/vj3LlzQv7ZZ59h4sSJkloREZEj+uyzzzB58mQha9u2LWJjY1GzZk05pcihcaCTNKdOnUJAQADy8vIsWfXq1XH06FF07txZYjMiInIUJ06cQPfu3YVjlu7u7oiNjUWHDh0kNiNHxiMuJE2HDh2wZcsWIcvPz4der0dmZqakVkRE5CgyMzPLfA7UJ598wnFOUnGgk1Tjxo3D7NmzhSwxMRFTp04FH9whIqKqoigKJk+ejMTERCGfM2cOxo4dK6kV0T084kLSFRYWIigoCLGxsUK+fv16vPLKK5JaERGRmq1fvx4LFy4UsoCAAERERKBatWqSWhHdw4FONuHKlSvQ6XRIS0uzZE5OTjh06BB69+4tsRkREalNREQE+vfvD5PJZMnq1q0Lo9GIp59+WmIzont4xIVsgre3N3bu3AmNRmPJTCYTQkNDcevWLYnNiIhITW7duoXQ0FBhnGs0GuzcuZPjnGwGBzrZjKFDh+L1118Xsps3byIsLEy48igREVFFlJSUlHnHz4oVKzBkyBBJrYhK40Anm7JixQoMHjxYyH788UcsX75cUiMiIlKLZcuWISIiQsiGDBlS6s4hItl4Bp1sTmpqKnQ6Ha5duybk//znPzFy5EhJrYiIyJ7985//xKhRo4SsSZMmMBqNqFevnqRWRGXjQCebFBMTg969e6O4uNiS1a5dG/Hx8WjRooXEZkREZG8uXrwIX19f4RobLi4uiIyMREBAgMRmRGXjEReySd26dcN7770nZL9eUCI/P19SKyIisjf3uwDehg0bOM7JZnGgk82aM2cOwsLChOz48eN48cUXJTUiIiJ7M2fOHJw4cULIxowZU+oieUS2hEdcyKbl5OQgICAAZ86cEfKtW7di6tSpkloREZE92Lp1K6ZNmyZk7du3R0xMDDw8PCS1Ino4DnSyeWfOnIG/vz9yc3MtmZubG6Kjo9G1a1eJzYiIyFYlJCQgMDAQhYWFlszDwwPHjh1D27ZtJTYjejgecSGb165dO2zdulXICgoKoNfrkZGRIakVERHZqoyMDBgMBmGcA/fuUec4J3vAgU52ITQ0FHPnzhWypKQkTJo0CWazWVIrIiKyNWazGRMnTkRSUpKQz5s3DyEhIZJaET0aHnEhu1FUVIQ+ffrg6NGjQr527VosXrxYUisiIrIla9euxdKlS4UsMDAQhw8fhqurq6RWRI+GA53syrVr16DT6ZCammrJtFotDhw4gL59+8orRkRE0h06dAgDBw4UHlmtV68eEhIS0LhxY4nNiB4Nj7iQXWnSpAl27doFjUZjycxmM8LCwnDz5k2JzYiISKbk5GSEhYUJ41yj0WD37t0c52R3ONDJ7gwaNAgrV64UspSUFISGhgpXHiUiIsdQXFyM0NBQ3L59W8hXrVqFgQMHSmpFVHEc6GSXli1bhmHDhgnZTz/9hNdee01SIyIikmXp0qWIjIwUsuHDh/N7AtktnkEnu3X37l3odDpcvXpVyL/88ks8//zzkloREZE1ffnllzAYDELm7e0No9EILy8vSa2IHg8HOtm1Y8eOoVevXigqKrJktWrVQlxcHFq1aiWxGRERVbXz58/Dz88P2dnZlszV1RVHjhyBn5+fxGZEj4dHXMiu+fv74/333xeyrKws6PV65OXlSWpFRERVLS8vDwaDQRjnAPDBBx9wnJPd40AnuzdjxgyMHz9eyE6ePInZs2eDDxAREamPoiiYNWsWTp48KeQTJkzACy+8IKkVUeXhERdShdzcXHTr1g2nTp0S8i1btmD69OmSWhERUVXYsmULZsyYIWQdO3ZETEwM3N3dJbUiqjwc6KQa586dg7+/f6mziFFRUfD19ZXYjIiIKktcXBx69uwpPPeoZs2aiIuLQ+vWrSU2I6o8POJCqtGmTRts27ZNyIqKimAwGJCWliapFRERVZa0tDQYDAZhnAPAp59+ynFOqsKBTqpiMBiwYMECIbt8+TImTpwoXF2OiIjsi9lsxoQJE3DlyhUhf/nll6HX6yW1IqoaHOikOuvWrUPPnj2FbO/evVi7dq2kRkRE9LjWrFmDb7/9Vsh69erFf9tJlXgGnVTpxo0b0Ol0wmWftVotvv/+ewwYMEBiMyIielT79+/HkCFDhEdCn3jiCSQkJOCpp56S2IyoanCgk2odPHgQgwYNEv5Br1+/PhISEtCoUSOJzYiIqLyuX78OHx8fpKamWjKtVov9+/ejX79+EpsRVR0ecSHV6t+/P9566y0hu3PnDkJCQlBcXCypFRERlVdRURFCQkKEcQ4Aq1ev5jgnVeNAJ1VbsmQJRowYIWRRUVFYtGiRpEZERFReixYtQnR0tJA9++yzWLx4saRGRNbBIy6keunp6fD19cWlS5eEfM+ePQgODpbUioiIHmTPnj0IDQ0VsmbNmiE+Ph516tSR1IrIOjjQySEYjUb06NEDhYWFlszDwwNxcXFo06aNxGZERPR7Z8+ehb+/P3JycixZtWrVEB0dDR8fH4nNiKyDR1zIIeh0OmzcuFHIcnJyoNfrkZubK6kVERH93q//Nv92nAPAhx9+yHFODoMDnRzGtGnTMGnSJCE7deoUZsyYAT6QREQkn6IomDFjBk6fPi3kkydPRnh4uKRWRNbHIy7kUPLy8hAYGIiff/5ZyDdt2oRZs2ZJakVERMC9f4vnzJkjZF26dEFUVBTc3d0ltSKyPg50cjiJiYnw9fVFVlaWJXN1dUVkZCT8/f0lNiMiclyxsbHo1auX8DK4tWrVQnx8PFq2bCmxGZH18YgLOZyWLVti+/btQlZUVASDwYC7d+/KKUVE5MBSU1MRHBxc6hoVn332Gcc5OSQOdHJIo0ePxquvvipkV69exfjx44UrjxIRUdUymUwYP348rl69KuQLFy7EqFGjJLUikosDnRzWmjVr0Lt3byHbt28fVq9eLakREZHjWb16Nb777jshCwoKwh//+EdJjYjk4xl0cmg3b96Ej48PUlJSLJlGo8G+ffswePBgic2IiNTvu+++w7Bhw4RX0mrQoAGMRiMaNmwosRmRXBzo5PB+/PFHDBgwACaTyZLVrVsXCQkJaNKkicRmRETqdfXqVeh0OuG5P05OTjh48CCCgoIkNiOSj0dcyOH16dOn1EOpd+/eRXBwMIqKiiS1IiJSr6KiIoSEhJR6Yv6aNWs4zonAgU4E4N6TkZ577jkhi4mJKfVEUiIienyvvPIKYmJihGzUqFH8N5fov3jEhei/MjIy4Ofnh4sXLwr57t27ERYWJqkVEZG67Nq1C+PGjROyFi1aIC4uDp6enpJaEdkWDnSi3zh+/DgCAwNRUFBgyWrUqIHY2Fi0b99eYjMiIvt3+vRp+Pv7Iy8vz5K5ubnh6NGj6NKli8RmRLaFR1yIfqNr167YtGmTkOXm5sJgMCAnJ0dSKyIi+5ednQ29Xi+McwDYvHkzxznR73CgE/3OlClTEB4eLmRnzpzB9OnTwQeciIgenaIomD59Os6ePSvk06ZNw+TJk+WUIrJhPOJCVIb8/Hz06NEDx48fF/KNGzfixRdflNSKiMg+bdy4EfPmzRMyHx8fREVFwc3NTVIrItvFgU50H0lJSdDpdMjMzLRkLi4uiIiIQPfu3SU2IyKyH9HR0QgKCkJJSYkl8/T0RHx8PJo3by6xGZHt4hEXovto3rw5Pv/8cyErLi5GcHAw7ty5I6kVEZH9uHPnDkJCQoRxDgCff/45xznRA3CgEz3AyJEjsWTJEiG7fv06xo0bJ1x5lIiIRCaTCWPHjsX169eFfOnSpXj22WcltSKyDxzoRA/x1ltvoV+/fkL2ww8/YNWqVZIaERHZvpUrV2L//v1C1q9fP/7bSVQOPINOVA4pKSnw8fHBzZs3LZlGo8HevXsxbNgwic2IiGzPt99+i2eeeUbIGjZsiISEBDz55JOSWhHZDw50onL66aef0K9fP+Foi5eXF4xGI7y9vSU2IyKyHZcvX4ZOp0N6erolc3JywuHDh9GrVy+JzYjsB4+4EJVT79698fbbbwtZWloaDAYDCgsLJbUiIrIdhYWFCA4OFsY5ALzzzjsc50SPgAOd6BEsWLAAer1eyOLi4rBgwQJJjYiIbMf8+fMRFxcnZAaDAfPnz5fUiMg+8YgL0SPKysqCn58fLly4IOQ7duzAuHHjJLUiIpJrx44dmDBhgpC1bt0ax44dQ61atSS1IrJPHOhEFXDy5El069YN+fn5lszd3R0xMTHo2LGjxGZERNZX1r+J1atXR0xMDDp16iSxGZF94hEXogro1KkTPvroIyHLy8uDXq9HVlaWpFZERNaXlZUFvV4vjHMA+PjjjznOiSqIA52ogiZOnIgZM2YI2fnz5xEeHg4+MEVEjkBRFEydOrXUkb+ZM2eWOu5CROXHIy5Ej6GgoAC9evVCfHy8kG/YsIFPiiIi1duwYQNefvllIfPz80NkZCSqVasmqRWR/eNAJ3pMZb3mr7OzMw4fPoyePXtKbEZEVHWOHDmCvn37oqSkxJLVqVMHRqMRTZs2lVeMSAV4xIXoMTVt2hR//etfhaykpAQhISG4ffu2pFZERFUnJSUFISEhwjgH7r2SC8c50ePjQCeqBM888wyWLVsmZMnJyRgzZoxw5VEiIntXUlKCMWPGIDk5WciXL1+O4cOHS2pFpC4c6ESVZOXKlRgwYICQHTx4ECtWrJDUiIio8q1YsQKHDh0SsoEDB+LNN9+UU4hIhXgGnagS3blzBz4+Prhx44aQf/PNNxgxYoSkVkREleObb77ByJEjhaxx48YwGo2oX7++pFZE6sOBTlTJoqOjERQUJJzN9PT0hNFoRLNmzSQ2IyKquKSkJPj6+iIjI8OSOTs7IyIiAoGBgRKbEakPj7gQVbLAwECsX79eyDIyMmAwGFBQUCCpFRFRxRUUFMBgMAjjHADeffddjnOiKsCBTlQF5s2bh5CQECEzGo2YN2+epEZERBU3d+5cJCQkCFlISAjmzp0rqRGRuvGIC1EVyc7Ohr+/P86dOyfk27dvx6RJkyS1IiJ6NNu3b8eUKVOErG3btoiNjUXNmjUltSJSNw50oip06tQpBAQEIC8vz5K5ubkhJiYGnTt3ltiMiOjhTpw4ge7duwvH82rUqIHY2Fi0b99eYjMideMRF6Iq1KFDB3zyySdCVlBQAL1ej8zMTEmtiIgeLjMzs8znzmzZsoXjnKiKcaATVbGxY8di9uzZQpaYmIgpU6aAD2ARkS1SFAWTJ09GYmKikM+ZMwdjx46V1IrIcfCIC5EVFBYWIigoCLGxsUL+zjvv4NVXX5XUioiobO+88w4WLVokZAEBAYiIiEC1atUktSJyHBzoRFZy5coV6HQ6pKWlWTInJyccPHgQQUFBEpsREf1PREQE+vfvD5PJZMnq1q0Lo9GIp59+WmIzIsfBIy5EVuLt7Y2dO3dCo9FYMpPJhNDQUNy6dUtiMyKie27evInQ0FBhnGs0GuzcuZPjnMiKONCJrGjo0KFYsWKFkN26dQthYWHClUeJiKytpKQEYWFhpe4weOONNzBkyBBJrYgcEwc6kZW9/vrrGDx4sJD9+OOPWL58uaRGRETAsmXLEBERIWRDhgzB66+/LqkRkePiGXQiCVJTU6HT6XDt2jUh//rrr/Hcc89JakVEjurrr7/G6NGjhaxJkyYwGo2oV6+epFZEjosDnUiSmJgY9O7dG8XFxZasdu3aiI+PR4sWLSQ2IyJHkpiYCF9fX2RlZVkyFxcXREZGIiAgQGIzIsfFIy5EknTr1g0bNmwQsszMTOj1euTn50tqRUSOJD8/HwaDQRjnAPCnP/2J45xIIg50Iolmz56NMWPGCNmJEyfw4osvSmpERI5kzpw5OHHihJCNHTsWs2bNktSIiAAecSGSLicnBwEBAThz5oyQb926FVOnTpXUiojUbuvWrZg2bZqQtW/fHrGxsahRo4akVkQEcKAT2YQzZ87A398fubm5lszNzQ1RUVHw8fGR2IyI1CghIQGBgYEoLCy0ZB4eHjh27Bjatm0rsRkRATziQmQT2rVrh61btwpZQUEBDAYDMjIyJLUiIjVKT0+HXq8Xxjlw7x51jnMi28CBTmQjQkNDMW/ePCFLSkrCpEmTYDabJbUiIjUxm82YNGkSLl26JOQvvfQSQkJCJLUiot/jERciG1JUVIS+ffsiOjpayNeuXYvFixdLakVEarF27VosXbpUyAIDA3H48GG4urpKakVEv8eBTmRjrl+/Dh8fH6SmployrVaLAwcOoG/fvvKKEZFdO3ToEAYOHCg8IlevXj0kJCSgcePGEpsR0e/xiAuRjWncuDF27doFjUZjycxmM8LCwpCcnCyxGRHZq+TkZISFhQnjXKPRYPfu3RznRDaIA53IBg0aNAgrV64UspSUFISGhgpXHiUiepji4mKEhITg9u3bQr5q1SoMHDhQUisiehAOdCIbtWzZMgwbNkzIIiMjS50fJSJ6kCVLluDIkSNCNnz4cLz22muSGhHRw/AMOpENu3v3Lnx9fXHlyhUh//LLL/H8889LakVE9uLLL7+EwWAQMm9vbxiNRnh5eUlqRUQPw4FOZOPi4uLQs2dPFBUVWbJatWohLi4OrVq1ktiMiGzZ+fPn4efnh+zsbEvm6uqKI0eOwM/PT2IzInoYHnEhsnF+fn54//33hSwrKwt6vR55eXmSWhGRLcvLy4PBYBDGOQB88MEHHOdEdoADncgOzJgxA+PHjxeykydPYtasWeCDYET0W4qiYObMmTh58qSQT5gwAS+88IKkVkT0KHjEhchO5Obmolu3bjh16pSQf/zxx/ymS0QWH3/8MWbOnClkHTt2RExMDNzd3SW1IqJHwYFOZEfOnTsHf3//UmdKo6Ki4OvrK7EZEdmCsp6zUrNmTcTFxaF169YSmxHRo+ARFyI70qZNG2zbtk3IioqKYDAYkJaWJqkVEdmCtLQ0GAwGYZwDwPbt2znOiewMBzqRnTEYDFiwYIGQXb58GRMnThSuEkhEjsNsNmPChAmlXpL15Zdf5kuyEtkhDnQiO7Ru3Tr07NlTyPbu3Yu1a9dKakREMq1ZswbffvutkPXq1Yv/JhDZKZ5BJ7JTN27cgE6nEy7frdVq8f3332PAgAESmxGRNe3fvx+DBw8WXtHpiSeeQEJCAp566imJzYioojjQiezYwYMHMWjQIOFoS/369ZGQkIBGjRpJbEZE1nD9+nX4+PggNTXVkmm1Wuzfvx/9+vWT2IyIHgePuBDZsf79+2P16tVCdufOHYSEhKC4uFhSKyKyhqKiIoSEhAjjHAD+8Ic/cJwT2TkOdCI7t3jxYowYMULIoqKisGjRIkmNiMgaFi1ahOjoaCF79tln+XefSAV4xIVIBdLT0+Hr64tLly4J+Z49exAcHCypFRFVlT179iA0NFTImjVrhvj4eNSpU0dSKyKqLBzoRCphNBrRo0cPFBYWWjIPDw8cO3YMbdu2ldiMiCrT2bNn4e/vj5ycHEtWrVo1REdHw8fHR2IzIqosPOJCpBI6nQ4bN24UspycHBgMBuTm5kpqRUSVKScnB3q9XhjnAPDhhx9ynBOpCAc6kYpMmzYNkydPFrJTp05hxowZ4INlRPZNURTMmDEDp0+fFvIpU6YgPDxcUisiqgo84kKkMnl5eQgMDMTPP/8s5Js2bcKsWbMktSKix7Vp0ybMmTNHyLp06YLo6GhUr15dUisiqgoc6EQqlJiYCF9fX2RlZVkyV1dXREZGwt/fX2IzIqqI2NhY9OrVS3j51Fq1aiE+Ph4tW7aU2IyIqgKPuBCpUMuWLbF9+3YhKyoqgsFgwN27d+WUIqIKSU1NRXBwcKlrG3z22Wcc50QqxYFOpFKjR4/Gq6++KmRXr17F+PHjhSuPEpHtMplMGD9+PK5evSrkCxcuxKhRoyS1IqKqxoFOpGJr1qxBUFCQkO3bt6/U1UeJyDatXr0a3333nZAFBQXhj3/8o6RGRGQNPINOpHI3b96Ej48PUlJSLJlGo8G+ffswePBgic2I6EG+++47DBs2THgFpgYNGiAhIQENGjSQ2IyIqhoHOpED+PHHHzFgwACYTCZLVrduXSQkJKBJkyYSmxFRWa5evQqdTic8Z8TJyQkHDx4s9agYEakPj7gQOYA+ffqUekj87t27CA4ORlFRkaRWRFSWwsJCBAcHl3pCd1lH1ohInTjQiRzEwoUL8dxzzwlZTEwMXnnlFUmNiKgsr7zyCmJjY4Vs1KhRpZ70TUTqxSMuRA4kIyMDfn5+uHjxopDv3r0bYWFhkloR0a927dqFcePGCVnLli0RFxeH2rVrS2pFRNbGgU7kYE6cOIHu3bujoKDAktWoUQOxsbFo3769xGZEju306dPw9/dHXl6eJXNzc8PRo0fRpUsXic2IyNp4xIXIwXTp0gWbNm0SstzcXBgMBuTk5EhqReTYsrOzodfrhXEOAJs3b+Y4J3JAHOhEDmjKlCkIDw8XsjNnzmD69Ongg2pE1qUoCqZPn46zZ88K+bRp0zB58mQ5pYhIKh5xIXJQ+fn56NGjB44fPy7kH3zwAebOnSupFZHj+eCDD/DSSy8JmY+PD6KiouDm5iapFRHJxIFO5MCSkpKg0+mQmZlpyVxcXBAREYHu3btLbEbkGKKjoxEUFISSkhJL5unpifj4eDRv3lxiMyKSiUdciBxY8+bN8fnnnwtZcXExgoODcefOHUmtiBzDnTt3EBISIoxzAPj88885zokcHAc6kYMbOXIklixZImTXr1/HuHHjhCuPElHlMZlMGDt2LK5fvy7kS5cuxbPPPiupFRHZCg50IsJbb72Ffv36CdkPP/yAVatWSWpEpG4rV67E/v37haxfv378O0dEAHgGnYj+KyUlBT4+Prh586Yl02g02Lt3L4YNGyaxGZG6fPvtt3jmmWeErGHDhkhISMCTTz4pqRUR2RIOdCKy+Omnn9CvXz/haIuXlxeMRiO8vb0lNiNSh8uXL0On0yE9Pd2SOTk54fDhw+jVq5fEZkRkS3jEhYgsevfujbffflvI0tLSYDAYUFhYKKkVkToUFhYiODhYGOcA8M4773CcE5GAA52IBAsWLIBerxeyuLg4LFiwQFIjInWYP38+4uLihMxgMGD+/PmSGhGRreIRFyIqJSsrC35+frhw4YKQ79ixA+PGjZPUish+7dixAxMmTBCy1q1b49ixY6hVq5akVkRkqzjQiahMJ0+eRLdu3ZCfn2/J3N3dERMTg44dO0psRmRfyvq7VL16dcTExKBTp04SmxGRreIRFyIqU6dOnfDRRx8JWV5eHvR6PbKysiS1IrIvWVlZ0Ov1wjgHgI8//pjjnIjuiwOdiO5r4sSJmDFjhpCdP38e4eHh4INvRA+mKAqmTp1a6qjYzJkzSx13ISL6LR5xIaIHKigoQK9evRAfHy/kGzZs4JPbiB5gw4YNePnll4XMz88PkZGRqFatmqRWRGQPONCJ6KHKeu1mZ2dnHD58GD179pTYjMg2HTlyBH379kVJSYklq1OnDoxGI5o2bSqvGBHZBR5xIaKHatq0Kf76178KWUlJCUJCQnD79m1JrYhsU0pKCkJCQoRxDtx7JReOcyIqDw50IiqXZ555BsuWLROy5ORkjBkzRrjyKJEjKykpwZgxY5CcnCzky5cvx/DhwyW1IiJ7w4FOROW2cuVKDBgwQMgOHjyIFStWSGpEZFtWrFiBQ4cOCdnAgQPx5ptvyilERHaJZ9CJ6JHcvn0bOp0ON27cEPJvvvkGI0aMkNSKSL5vvvkGI0eOFLLGjRvDaDSifv36kloRkT3iQCeiRxYdHY2goCDhjK2npyeMRiOaNWsmsRmRHElJSfD19UVGRoYlc3Z2RkREBAIDAyU2IyJ7xCMuRPTIAgMDsX79eiHLyMiAwWBAQUGBpFZEchQUFCA4OFgY5wDw7rvvcpwTUYVwoBNRhcybNw8hISFCZjQa8dJLL0lqRCTHvHnzYDQahSwkJARz586V1IiI7B2PuBBRhWVnZ8Pf3x/nzp0T8u3bt2PSpEmSWhFZz/bt2zFlyhQha9u2LWJjY1GzZk1JrYjI3nGgE9FjOXXqFAICApCXl2fJ3NzcEBMTg86dO0tsRlS151ez+wAAHqxJREFUTpw4ge7duwvHumrUqIHY2Fi0b99eYjMisnc84kJEj6VDhw745JNPhKygoAB6vR6ZmZmSWhFVrczMzDKfc/HJJ59wnBPRY+NAJ6LHNnbsWMyePVvIEhMTMXXqVPBBOlIbRVEwefJkJCYmCvmcOXMwZswYSa2ISE14xIWIKkVhYSGCgoIQGxsr5OvXr8crr7wiqRVR5Vu/fj0WLlwoZAEBAYiIiEC1atUktSIiNeFAJ6JKc+XKFeh0OqSlpVkyJycnHDx4EEFBQRKbEVWOiIgI9O/fHyaTyZLVrVsXRqMRTz/9tMRmRKQmPOJCRJXG29sbO3fuhEajsWQmkwmhoaG4deuWxGZEj+/mzZsIDQ0VxrlGo8HOnTs5zomoUnGgE1GlGjp0KFasWCFkt27dQlhYmHDlUSJ7UlJSgrCwsFI/aL7xxhsYMmSIpFZEpFYc6ERU6V5//XUMHjxYyH788UcsX75cUiOix7Ns2TJEREQI2ZAhQ/D6669LakREasYz6ERUJVJTU6HT6XDt2jUh//rrr/Hcc89JakX06L7++muMHj1ayJo0aQKj0Yh69epJakVEasaBTkRVJiYmBr1790ZxcbElq127NuLj49GiRQuJzYjKJzExEb6+vsjKyrJkLi4uiIyMREBAgMRmRKRmPOJCRFWmW7dueO+994QsMzMTer0e+fn5kloRlU9+fj4MBoMwzgFgw4YNHOdEVKU40ImoSpV18ZYTJ07gxRdflNSIqHzmzJmDEydOCFlZF+UiIqpsPOJCRFUuJycHAQEBOHPmjJBv3boVU6dOldSK6P62bt2KadOmCVn79u0RGxuLGjVqSGpFRI6CA52IrOLMmTPw9/dHbm6uJXNzc0N0dDS6du0qsRmRKCEhAYGBgSgsLLRkHh4eOHbsGNq2bSuxGRE5Ch5xISKraNeuHbZu3SpkBQUF0Ov1yMjIkNSKSJSeng6DwSCMc+DePeoc50RkLRzoRGQ1oaGhmDt3rpAlJSVh0qRJMJvNkloR3WM2mzFp0iQkJSUJ+bx58xASEiKpFRE5Ih5xISKrKioqQp8+fXD06FEhX7t2LRYvXiypFdG9r8GlS5cKWWBgIA4fPgxXV1dJrYjIEXGgE5HVXb9+HT4+PkhNTbVkWq0WBw4cQN++feUVI4d16NAhDBw4UHgkp169ekhISEDjxo0lNiMiR8QjLkRkdY0bN8auXbug0WgsmdlsRlhYGG7evCmxGTmi5ORkhIWFCeNco9Fg9+7dHOdEJAUHOhFJMWjQIKxcuVLIUlJSEBISIlx5lKgqFRcXIzQ0FLdv3xbyVatWYeDAgZJaEZGj40AnImmWLVuGYcOGCVlkZGSpc8BEVWXJkiWIjIwUsuHDh+O1116T1IiI/r+9O4+Oos73Pv7pdBJIgAhRwIU9UZDlShYSEYwg7iCjZjGyjYISeFBEZtTx0fEZZ+Ig4zA+ChfEbZTB5TIyc89l8EHJgxuQCSRBiISIJCCLCGHNvnSn7h9cM5YhQpLururO+3WOR86X7q5POpzKJ9W/qgJr0AFY7Pjx44qLi9M333xjmq9evVp33XWXRanQHqxevVopKSmmWd++fZWfn6/IyEiLUgEABR2ADeTm5mrUqFGqq6trnEVERCg3N1eXX365hckQqHbv3q34+HiVl5c3zkJDQ7Vp0ybFx8dbmAwAWOICwAbi4+P14osvmmZlZWVKTk5WVVWVRakQqKqqqpSSkmIq55L00ksvUc4B2AIFHYAtZGRkaMqUKaZZQUGBZs+eLT7og6cYhqFZs2apoKDANJ86dapmzpxpUSoAMGOJCwDbqKysVGJionbu3GmaL1++nPIEj1i+fLlmzZplmg0dOlQ5OTkKDw+3KBUAmFHQAdjKV199pREjRjRZG7x582bFxcVZmAz+7mznOnTp0kW5ubm64oorLEwGAGYscQFgKwMHDtQbb7xhmtXV1SklJUUnTpywKBX83YkTJ5SSkmIq55L05z//mXIOwHYo6ABsJyUlRY888ohptm/fPk2bNs10t0fgfDQ0NGjq1KlNLuU5f/58JScnW5QKAJpHQQdgSwsXLtSoUaNMs7Vr1+q5556zKBH81YIFC/TBBx+YZqNHj+bfEgDbYg06ANs6dOiQYmNjTbdhDwoK0kcffaRx48ZZmAz+IisrSzfffLPpk5cePXpo27ZtuvTSSy1MBgDNo6ADsLUNGzboxhtvNBWs7t27a9u2bbrsssssTAa7O3jwoGJiYnTs2LHGWVBQkLKysjR27FgLkwHAT2OJCwBbu/7665WZmWmalZaWKi0tTfX19Ralgt3V1dUpLS3NVM4l6dlnn6WcA7A9CjoA23v88cc1YcIE02zz5s167LHHLEoEu3vssceUnZ1tmt1+++38mwHgF1jiAsAvnDx5UnFxcdq7d69pvmrVKqWmplqUCna0atUq3X333aZZ//79lZeXp27dulmUCgDOHwUdgN/Iz8/XNddco9ra2sZZ586dlZubq4EDB1qYDHZRVFSkESNGqKKionHWoUMHZWdnKyYmxsJkAHD+WOICwG/ExsZq8eLFpllFRYWSk5NVWVlpUSrYxff/Fn5YziVpyZIllHMAfoWCDsCv3H///br33ntNs507dyojI0N8INh+GYahjIwMFRYWmub33XefZsyYYVEqAGgdlrgA8DtVVVUaOXKkduzYYZovXbpUs2fPtigVrLR06VLNmTPHNLvqqquUnZ2tsLAwi1IBQOtQ0AH4pT179iguLk5lZWWNs9DQUG3cuFEjRoywMBl8bcuWLRo9erTpspsRERHKy8tTdHS0hckAoHVY4gLAL0VHR+vNN980zerq6pSSkqLjx49bEwo+d+zYMaWmpja5Jv5bb71FOQfgtyjoAPzWnXfeqV/+8pem2f79+zVlyhTTnUcRmNxut6ZMmaL9+/eb5o8++qjuuOMOi1IBQNtR0AH4tQULFujaa681zdatW9fk7qMIPJmZmfrwww9Ns6SkJP3+97+3KBEAeAZr0AH4vcOHDysmJkZHjhxpnDkcDq1bt0433XSThcngLR9++KFuvfVW05V7Lr74YuXn5+uSSy6xMBkAtB0FHUBA+PTTTzVu3Di53e7G2YUXXqht27apd+/eFiaDp+3fv1+xsbGmcw2cTqc2bNigpKQkC5MBgGewxAVAQLjuuuuaLG04fvy4UlNTVVdXZ1EqeFpdXZ3S0tKanAi8YMECyjmAgEFBBxAwHn30Uf3sZz8zzXJycpqcSAr/9Ytf/EI5OTmm2R133MH3GEBAYYkLgIBy6tQpxcfHq7i42DR/9913lZ6eblEqeMI777yjyZMnm2ZRUVHKzc1V165dLUoFAJ5HQQcQcL744guNHDlSNTU1jbNOnTppy5YtGjx4sIXJ0FqFhYUaMWKEqqqqGmcdO3bUP//5T1111VUWJgMAz2OJC4CAM3z4cC1dutQ0q6ysVEpKiioqKixKhdYqLy9XcnKyqZxL0rJlyyjnAAISBR1AQLrvvvs0Y8YM02zXrl164IEHxAeH/sMwDD3wwAMqKioyze+//37de++91oQCAC9jiQuAgFVdXa1rrrlGX3zxhWm+ePFiPfjggxalQkssXrxYc+fONc1iYmK0efNmdezY0aJUAOBdFHQAAa2kpESxsbE6ffp04ywkJESfffaZrr76aguT4Vyys7OVlJQkl8vVOOvatavy8vI0YMAAC5MBgHexxAVAQBswYIBWrFhhmtXX1ys1NVWlpaUWpcK5lJaWKi0tzVTOJWnFihWUcwABj4IOIOBNnDhRv/rVr0yzgwcPavLkyaY7j8Ie3G63Jk2apIMHD5rmTzzxhG6//XaLUgGA71DQAbQLv/vd7zR27FjTbP369frtb39rUSI055lnnlFWVpZpNnbsWL5XANoN1qADaDeOHDmimJgYHT58uHHmcDi0du1a3XrrrRYmw/c++OADjR8/3jS75JJLtG3bNvXs2dOiVADgWxR0AO3K559/rrFjx5qWtkRGRio/P199+/a1MBn27dun2NhYnTx5snHmdDr1ySefaPTo0RYmAwDfYokLgHbl2muv1cKFC02zEydOKCUlRbW1tRalQm1trVJTU03lXJL+8Ic/UM4BtDsUdADtzvz583XXXXeZZrm5uXrkkUcsSoR58+YpNzfXNEtOTuZ7AqBdYokLgHaprKxM8fHx+vrrr03zlStXavLkyRalap9WrlypqVOnmmZXXHGFtm7dqoiICItSAYB1KOgA2q2CggIlJiaqurq6cRYeHq4tW7ZoyJAhFiZrP7788kslJCSYvgdhYWHKycnRsGHDLEwGANZhiQuAdmvYsGF6+eWXTbOqqirdddddKisrsyhV+1FWVqbk5GRTOZek5cuXU84BtGsUdADt2rRp0zRz5kzTbPfu3ZoxY4b4gNF7DMPQ9OnTtXv3btM8IyOjyXIXAGhvWOICoN2rqanR6NGjlZeXZ5q/8MILmjdvnkWpAtsLL7yg+fPnm2bx8fH6/PPP1bFjR4tSAYA9UNABQGe/BndwcLA++eQTjRo1ysJkgWfTpk0aM2aMXC5X46xbt27Kz89Xv379rAsGADbBEhcAkNSvXz/95S9/Mc1cLpfS0tJ09OhRi1IFnqNHjyotLc1UzqUzV3KhnAPAGRR0APgf48eP15NPPmmaffvtt7rnnntMdx5F67hcLqWnp+vbb781zZ966inddtttFqUCAPuhoAPADzzzzDMaN26cabZhwwY9/fTTFiUKHE8//bQ+/vhj0+yGG27Qb37zG2sCAYBNsQYdAH7k6NGjio2N1aFDh0zzNWvWaMKECRal8m9r1qzRxIkTTbNevXopPz9f3bt3tygVANgTBR0AziI7O1tJSUmmtdJdu3ZVfn6++vfvb2Ey/1NSUqK4uDidOnWqcRYcHKzPPvtMI0eOtDAZANgTS1wA4CxGjhypP/7xj6bZqVOnlJKSopqamp98rmEYOnDggOrq6rwZ0VJ1dXU6cODAOa8VX1NTo9TUVFM5l6RFixZRzgGgGRR0AGjG3LlzlZaWZprl5+fr4YcfbvY5BQUFGj58uPr06aPBgwdr3759Xk7pe3v37tXgwYPVp08fDR8+XAUFBc0+du7cucrPzzfN0tLS9NBDD3k7JgD4LQo6ADTD4XDotdde08CBA03zV155RW+99VaTx69YsUKJiYnasWOHJKm4uFiLFi3ySVZfWrRokYqLiyVJO3bsUGJiYpNLVErSm2++qVdffdU0GzRokF577TU5HA6fZAUAf0RBB4Cf0KVLF61evVrh4eGm+axZsxqLeE1NjTIyMvTzn/9c1dXVpsdt2rTJZ1l95cdfU3V1taZNm6ZZs2Y1Lv/Zvn27Zs+ebXpceHi43n//fXXp0sVnWQHAH3GSKACch7fffltTpkwxzaKjo7V69WpNnz5deXl5Z31eSEiIysvL1aFDB1/E9Lra2lp17ty5yY2GvhcXF6c33nhDycnJ2rNnj+nv3n77bU2aNMkXMQHAr1HQAeA8zZkzR0uXLjXNQkJCVF9f/5PPy8vLU2xsrDej+UxeXp7i4+N/8jFne0/mzJmjJUuWeDMaAAQMlrgAwHn605/+pISEBNPsXOVcUpOTJP3Z+XwtP35PEhISAnItPgB4CwUdAM5Thw4dtGzZMgUHB7foee2toP9QSEiIXn755YBZ4gMAvkBBB4DzlJ2drYkTJza7/ro57bmg19fXa+LEicrOzvZSIgAIPBR0ADgHwzD04osvKikpSYcOHWrx87dv397iUm9H9fX12r59e4ufd/DgQSUlJemll146542NAAAUdAD4SeXl5UpPT9e8efNaXbJrampUVFTk4WS+V1RUpNra2lY91+Vy6eGHH1Z6errKy8s9nAwAAgsFHQCaUVpaqoSEBK1atarNrxUIy1w88TWsWrVKCQkJKi0t9UAiAAhMFHQAaEZmZqbHjnxT0P+lqKhImZmZHnktAAhEFHQAaEZlZaXHXouCbubJ9xYAAg03KgKAZhQVFenGG2/UwYMH2/xanTt31unTpxUU5J/HRRoaGhQREeGRYt2rVy+tX79egwYN8kAyAAg8/vmTAgB8YNCgQdq1a5eef/559ezZs02vVVFRoT179ngome99/fXXbS7nPXv21PPPP69du3ZRzgHgJ3AEHQDOQ3V1tV5//XUtXLiw1UfU3333XaWnp7cph7vB0MEyQ8Unzvy/1i3VuaV695m/D3FKoU6pg1PqFeFQVKRDvSIccgY52rTdd999V5MmTWrVc3v37q3HH39c06dPV1hYWJtyAEB70LLb4QFAOxUWFqYHH3xQM2fO1IoVK7RgwQKVlJS06DXy8/NbVNANw9D+02fKePFJQyUnDO09ZajO3bLsoU6pf1eHBkQ6FNXtTGnvc4FDDsf5l/bWrD+PiorSE088oalTpyo0NLTFzweA9ooj6ADQCi6XS++9956effbZ877Sy7hx45SVlXXOx1XVG/p0X4PW7WnQgdPe2UX3vsChW6KDdF2/IIWHnLuojxs3Ths2bDiv177yyiv15JNP6u6771ZwMMeBAKClKOgA0AYNDQ3629/+pszMzHPeZTMsLEyVlZXNHrn+5tSZUv7pvgbV+OjGox2Dpev6BemW6CD17Xr205IMw1B4eLhqamp+8rWGDx+up556SnfeeaffngwLAHZAQQcADzAMQ2vXrlVmZqZycnKafVxhYaGuvPJK0/O2HDK05iu3Ckut3R0P7u7Q7QOdSrjMvPxl586dGjp0aLPPS0xM1K9//WvddtttLVo2AwA4Ow5xAIAHOBwOTZgwQdnZ2crKytKYMWPO+rhTp0796881hv6wyaWFG12Wl3NJKiw1tHCjS89vculUzb/ynD59+qyPHzNmjLKyspSdna3x48dTzgHAQziCDgBesnHjRmVkZKiwsFCSNGDAgMZLLW7c36BX89yqqLMyYfO6hEoPxDk1qs+Z4zjR0dGNJ8UOHjxYy5cv1+jRo62MCAABi4IOAF6Wk5OjI0eOaMKECSqrc2h5rks5B/1j13t1L4dmxgcrItTQP/7xD1188cVKSEiwOhYABDQKOgD4SPaBBi3b6rLtUfPmdAmVZo0I1sjerIoEAF+goAOAD6zd7dbr+S28gLnNzIh1avwVTqtjAEDAo6ADgBcZhqG/7mzQe1/6dzn/XvpQp1KHBHFCKAB4EQUdALzonR0uvV/YYHUMj0oZHKRJ/8YNiADAW1hQCABe8p+73AFXziXp/cIG/WdRYHwiAAB2REEHAC9YX+zWiu2BW2JXfOHW+uLA/foAwEoUdADwsJKTDXolN/DL6yu5bu09GXifEACA1ViDDgAeVO829Nh6l7451bpda2XpPn00d6CHU/20QclP6cqUX7fquf26OrTwxmCFODlpFAA8hSPoAOBBqwvdrS7n/mjfKUOrCwP/0wIA8CUKOgB4SMnJBq0OwJNCz2V1YQNLXQDAg7hOFgB4QL3b0OIct9xtPHge5AzVBX3/zTOhfqTySIlcNRVN5uEX9WnT67oNaXGOWwtvdLDUBQA8gDXoAOAB7xW4tGqnfY8iH/9qsz7/3U0y3PWm+WUjU5Uwd6VHtpE2JEjpwzjuAwBtxRIXAGij8lpDfy+ybzmvPnFIOS+kNynnEX2GKTbjFY9t5+9FDSqv5ZgPALQVBR0A2mjD3gbV2/Q8SXd9rXL+dLdqTx8xzUM6R+rqX/xVwR3CPbaterf08V77/qICAP6Cgg4AbdBgGPpwj03buaQvXn9IJ4u3mmaOIKcS5q5Upx79Pb69D/e41cDKSQBoEwo6ALTB9u8Mfdf0vEtbKF7379r/6VtN5kPuyVSPYeO8ss3DFdKO7yjoANAWFHQAaIN1Nj16fqzwMxWsfKzJvNc1d+vyCfO9um27vicA4C8o6ADQSkcrDeUest/R4qrjB5Tz4iQZbpdpfkHfqxQz82Wvb3/rIUOllfZ7XwDAX1DQAaCV/n+JW3aroe66GuUsSlNdWalpHtr5QiX+YpVHTwptjqEz7w0AoHUo6ADQSjuP2q2eS9tena1Te/NNM0eQUyMeXqlO3fv5LIcd3xsA8BcUdABohQbD0N6T9iqhX6/9vzqw8Z0m8yGTfq8eQ6/3aZaSkwZXcwGAVqKgA0ArHC6Xql3nfpyvHP1yg3a+87+bzHuPukeXj5/n8zzVLum7cp9vFgACAgUdAFphzwn73JCnsnSftr44RUaDed33Bf2GK2bmMotS2es9AgB/QkEHgFYoscnyFldtlXIWpaqu4rhpHtrlIl09f5WcoWEWJZOKbfIeAYC/oaADQCsUn7BH+cxfPlOnv9lhmjmcwUp4+G2Fd+9rUaozSmzyHgGAv6GgA0ALGTY5QXT3f/1Rh7L/2mQ+dNICdR8yxveBfqTkpCGDE0UBoMUo6ADQQvUN1p8gemTHehX+x9NN5r2vnazo2+ZakKipapfkYhk6ALQYBR0AWqje4nvwVBwp1taXpjY5KbRr/1jF3L/UolRnV8f9igCgxSjoANBCtRaWTldNhXIWpaq+8qRpHhrRXYnz/0PO0I4WJTs7CjoAtBwFHQBayMplG3nL7lfZgZ2mWeNJoRf1sShV81jiAgAtR0EHgBYKsWjP+dXfn9O3W/7eZD5sykJ1H3ydBYnOLcRpdQIA8D8UdABooVALSud32/6fCv/6TJN5n6SpirrlQd8HOk9W/TIDAP6MXScAtJCvC3r54d3KXfJzyTCvF+k6IE7DZyzxbZgWsuKXGQDwdxR0AGih4CApPMQ326qvLj9zUmjVadO8Q0QPW54U+kOdQs68VwCAlmHXCQAt5HA4NKCbw+vbMQxDeUvvU/mhIvP2ncFKmPeOwi/s7fUMbTGgm0MOh/ffJwAINBR0AGiFqEjvF8+i1Zk6nLumyXzY1Od10ZXXen37beWL9wgAAhEFHQBaIcrLR9AP565R0d+ebTLvc900Rd38v7y6bU8ZQEEHgFahoANAK0RFem/3WX6oSLlLp0uGYZp3ixph+5NCfyi6Gz9iAKA12HsCQCtc3Nk7J4rWV53WPxelyFVdZpp3uKCnEh95T86QDp7fqBd0CpF6drY6BQD4p2CrAwCAP/r+RNEvjxrnfnALbH9zvioOf91kbhgNyn7+To9u64euf26rR1+PE0QBoPUo6ADQSkN7BOnLo26Pvmb1sf1nndeVlaqurNSj2/KmoT35gBYAWos9KAC00g1RQQriIHETToc0bgA/XgCgtdiDAkArRYY5lNiLhv5jib2CFBnG+wIAreUwDMOzCygBoB0pONKg//Oxy+oYtvLbscEscQGANmAPCgBtMLSHQ5dFWJ3CPnpFSEN6cPQcANqCgg4AbeBwOHRLtNPqGLZxS7STq7cAQBtR0AGgjcb0C1JHromljsHSdf34sQIAbcWeFADaqFOoQylDOIqeOsSpTqEcPQeAtqKgA4AH/GxgkC6PbL/l9PILHZo4kB8pAOAJ7E0BwAOcQQ49mBis4Ha4Vw0Jkh5KCJaTi8IDgEe0wx8lAOAdvS9wKH1Y+1vqkj7MqV4XUM4BwFMo6ADgQe1tqcvlkSxtAQBPY68KAB7Unpa6hARJDyWytAUAPK0d/AgBAN/qfYFDcxOdCuTa6pD0UCJLWwDAGyjoAOAFo/s6NTM+cNejZ8Q7Nbpv4H59AGAlCjoAeMnN0U5NvSrwSuzUq5y6ibunAoDXOAzDMKwOAQCB7B9fufXGNrfVMTxiRqxT46+gnAOAN1HQAcAHPt7r1pIct/x1hxvkkOYkODW2P+UcALyNgg4APvLF4QYt2eLSiWqrk7RMZJj0YEKwhl/CqkgA8AUKOgD4UGWdoTe2ufXx3garo5yX6/sH6b4YpzqFcrUWAPAVCjoAWCDv2wYt22rfo+mRYdLsEcGKu5Sj5gDgaxR0ALCIXY+mc9QcAKxFQQcAixUcadCar9zK+9aw7CRSh6S4Sx26faBTw3py1BwArERBBwCbOFph6KNit7JKGlRW65ttRnSQbhgQpJuinOrRmSPmAGAHFHQAsJl6t6HsAw1at6dBRce8s4sedJFDt0QHaWTvIIU4KeYAYCcUdACwsWNVhopPGCo+0aDik2f+3NKj6xEdpKhIh6K6ORQVGaSoSIcuCqeUA4BdUdABwI8YhqHjVVLxSUMHywzVugzVuaW6/7lRaajzzH8dgh3qFXGmlF8YLjkcFHIA8BcUdAAAAMBGOFUfAAAAsBEKOgAAAGAjFHQAAADARijoAAAAgI1Q0AEAAAAboaADAAAANkJBBwAAAGyEgg4AAADYCAUdAAAAsBEKOgAAAGAjFHQAAADARijoAAAAgI1Q0AEAAAAboaADAAAANkJBBwAAAGyEgg4AAADYCAUdAAAAsBEKOgAAAGAjFHQAAADARijoAAAAgI1Q0AEAAAAboaADAAAANkJBBwAAAGyEgg4AAADYCAUdAAAAsBEKOgAAAGAjFHQAAADARijoAAAAgI1Q0AEAAAAboaADAAAANvLf3BA1W2KD64kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 900x660 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(3, 2.2))\n",
    "G.plot(ax=ax, legend=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate our fake data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(20, 30)\n",
    "y = np.linspace(100, 300, 50)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "X = X.flatten()\n",
    "Y = Y.flatten()\n",
    "\n",
    "Z = 2*X+3*Y\n",
    "\n",
    "model = pd.DataFrame({'X': X, 'Y': Y, 'Z': Z})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And inspect it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>Z</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20.000000</td>\n",
       "      <td>100.0</td>\n",
       "      <td>340.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>20.204082</td>\n",
       "      <td>100.0</td>\n",
       "      <td>340.408163</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>20.408163</td>\n",
       "      <td>100.0</td>\n",
       "      <td>340.816327</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>20.612245</td>\n",
       "      <td>100.0</td>\n",
       "      <td>341.224490</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>20.816327</td>\n",
       "      <td>100.0</td>\n",
       "      <td>341.632653</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2495</th>\n",
       "      <td>29.183673</td>\n",
       "      <td>300.0</td>\n",
       "      <td>958.367347</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2496</th>\n",
       "      <td>29.387755</td>\n",
       "      <td>300.0</td>\n",
       "      <td>958.775510</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2497</th>\n",
       "      <td>29.591837</td>\n",
       "      <td>300.0</td>\n",
       "      <td>959.183673</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2498</th>\n",
       "      <td>29.795918</td>\n",
       "      <td>300.0</td>\n",
       "      <td>959.591837</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2499</th>\n",
       "      <td>30.000000</td>\n",
       "      <td>300.0</td>\n",
       "      <td>960.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2500 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "              X      Y           Z\n",
       "0     20.000000  100.0  340.000000\n",
       "1     20.204082  100.0  340.408163\n",
       "2     20.408163  100.0  340.816327\n",
       "3     20.612245  100.0  341.224490\n",
       "4     20.816327  100.0  341.632653\n",
       "...         ...    ...         ...\n",
       "2495  29.183673  300.0  958.367347\n",
       "2496  29.387755  300.0  958.775510\n",
       "2497  29.591837  300.0  959.183673\n",
       "2498  29.795918  300.0  959.591837\n",
       "2499  30.000000  300.0  960.000000\n",
       "\n",
       "[2500 rows x 3 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And to go the other way, we simply perform an OLS Linear Regression fit to estimate the dependency of Z as a function of X and Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols = smf.ols('Z ~ X + Y', model).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From where we easily extract the coefficients of our model for $f_Z$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>Z</td>        <th>  R-squared:         </th>  <td>   1.000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   1.000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>4.991e+30</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sat, 25 Jul 2020</td> <th>  Prob (F-statistic):</th>   <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>17:37:38</td>     <th>  Log-Likelihood:    </th>  <td>  62959.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  2500</td>      <th>  AIC:               </th> <td>-1.259e+05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  2497</td>      <th>  BIC:               </th> <td>-1.259e+05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>      <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>      <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>-1.012e-11</td> <td> 5.15e-13</td> <td>  -19.659</td> <td> 0.000</td> <td>-1.11e-11</td> <td>-9.11e-12</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X</th>         <td>    2.0000</td> <td>  1.9e-14</td> <td> 1.05e+14</td> <td> 0.000</td> <td>    2.000</td> <td>    2.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Y</th>         <td>    3.0000</td> <td>  9.5e-16</td> <td> 3.16e+15</td> <td> 0.000</td> <td>    3.000</td> <td>    3.000</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>1971.543</td> <th>  Durbin-Watson:     </th> <td>   0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td> 148.925</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td>-0.003</td>  <th>  Prob(JB):          </th> <td>4.58e-33</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 1.804</td>  <th>  Cond. No.          </th> <td>1.93e+03</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.93e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      Z   R-squared:                       1.000\n",
       "Model:                            OLS   Adj. R-squared:                  1.000\n",
       "Method:                 Least Squares   F-statistic:                 4.991e+30\n",
       "Date:                Sat, 25 Jul 2020   Prob (F-statistic):               0.00\n",
       "Time:                        17:37:38   Log-Likelihood:                 62959.\n",
       "No. Observations:                2500   AIC:                        -1.259e+05\n",
       "Df Residuals:                    2497   BIC:                        -1.259e+05\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept  -1.012e-11   5.15e-13    -19.659      0.000   -1.11e-11   -9.11e-12\n",
       "X              2.0000    1.9e-14   1.05e+14      0.000       2.000       2.000\n",
       "Y              3.0000    9.5e-16   3.16e+15      0.000       3.000       3.000\n",
       "==============================================================================\n",
       "Omnibus:                     1971.543   Durbin-Watson:                   0.000\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              148.925\n",
       "Skew:                          -0.003   Prob(JB):                     4.58e-33\n",
       "Kurtosis:                       1.804   Cond. No.                     1.93e+03\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.93e+03. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now build the model for SCM 1.5.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "G2 = CausalModel()\n",
    "G2.add_causation(r'$U_X$', 'X')\n",
    "G2.add_causation(r'$U_Y$', 'Y')\n",
    "G2.add_causation(r'$U_Z$', 'Z')\n",
    "G2.add_causation('X', 'Y')\n",
    "G2.add_causation('Y', 'Z')\n",
    "G2.pos = {r'$U_X$': (-1, 1), r'$U_Y$':(0,1), r'$U_Z$': (1,1), 'X': (-1, 0), 'Y':(0,0), 'Z':(1, 0)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAYAAACmKP9/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAuIwAALiMBeKU/dgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hb1eH/8c+VLNtx4hlnOstxBklIyA6jJGGkpEDZpQtKSxiBMkqh/bZQaBktdAC/FiiU2dIBHRQ6KGGHTSYBQhKSOE7iLDI8Y8e2xv39IWyieEmW7Ht19H49T5+iE+nqWDo6n3vOPfdey7ZtWwAAwCgepysAAAASj4AHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMlOZ0BZzSFLS1pcrWpkpblQds+UNSU1AKhKQ0j+TzSOleKTfT0sh8SyPyLGWkWU5XG+iUbdvaWy+VVtraVm2rKWirMSj5g+F/93mlDK+U7rU0JNdSSb6lwizJsmjfcL/GgK3Nn/bd1Q22moKSP/RZ353uDfff+b3CfffwPEvp3tRs2ykR8P5guEGUVtgqrQz/f3m1raAd/TY8ljQkx1JJQbhDLCkg9OEOe+tsbay0taki1NK+axpj20ZOhlra9sgCj0blWyrsTduGs5rD/OC+e1uNrVAMfbfXkobmtu67fSkQ+pZt2zF8VMllR62t5zcG9WpZSPubEr/9XmnS3GKP5o/yamiu+Y0F7uEP2nq7PKRFG0P6eG/3/ITHFlqaP8qjo4d6UqIzhHuUV9tatDGoxWUhHQgkfvt90qXjij06aZRXg7PNbdvGBXwwZGvFjnDjWLWr5/60w/tbmj/Kq5lDLKV5zG0wcNbu/bZeKA3qpU2hmEfpXZWTIZ040qPPl3jVvw9tG90jELK1dFu47169u+f67skDLc0f7dW0QZa8hvXdxgR8vd/W/9aH9EJpUHvrnatHfqY0r8SjU8Z4lZ1hVmOBcz78JKT/fBzUih22nPrBWpKmDbb0xbFeTRzA+lwkRm2jrWfXB/ViaUiVDc7VozBL+nyJVyeP8SjLZ0bfbUTAr9oZ0m+XBRwN9kPlZUoLp6dp5hA6QnRdXZOtR98LH2Zyk+OKPbpwile9083oCOGMpdtCemB5QFUOBvuhCrOky2ekafKg5O+7kzrg6/22fv9eeLrSrWYP92jBVEbziN2KHSHdvyygigNO16RtBb2ky2akadrg5O8I0bNqG209vDKoN7a4t++eN9KjC6Z4k3o0n7QB78ZRe3vyMqWFM9I0s4iOEJ1z66i9PYzmEQs3jtrbk+yj+aQL+IaArcdWBvWii0ft7Zkz3KOLptERon2rdoZ071L3jtrbU9BLumJm8naE6H51TbYeWhHU6y4etbdn3kiPLpzqTbrTopMq4OuabN32eqDbTgvqCcX5lm6ak6bczORqKOh+r2wK6r6lQccW0cXLknTFLK+OK/Y6XRW4TFWDrVtfC6isMllbt3RYoaUbZqcl1QAtaQK+qsHWrYsDKqtKiup2qChH+vEcHxcSQYv/fhzUo+8Fna5GQiyY6tUpYwh5hO2ts3XzYr+21zpdk/gV51m6cW6a8pJkgJYUAV/VYOvGl81oIM36ZUk/PYGQh/T02qD++L4Z4d7s/CO8OnMcIZ/q9tTZ+tHLfu1JgrVS0SrKlm49wZcUIe/6A2Z1TbZuWRwwKtwlaU+9dPNiv6oaXL9/hW70/Ebzwl2S/vh+UM9vNO/vQvSqGsIjd5PCXZK210q3Lg6orsn9fberA74hYOunrwe02YBp+bZsr5VufS05GgoS740tQT243NwQfHB5UG9uMffvQ/vqmsKHVHcYNjBrVlYVzqaGgLv7blcH/GMrg1qXxAvqolFWaesBgzt5tK282tY9S5J3QV00bEm/WRJUebXJfyXacv+yoBHrpTqybq+tx1y+bsa1d5NbtTPU7qlwW9/4s1b89sKIshHHX6gpF98f9/tuevF3ev/RqyLKRp/6XR3+9dvj3nZ73toa0tFDQzpqqKv3t5AgwZCte5YEFGjnbCGT2ncgJN27JKCfnZhm3HW+0ba3y0N6uzw1+u4XS0M6amhIkwe6s+92Za3qmmz9dln7txCqKnuvVVneiCkJee82t12cmG135HfLA6rmeHxK+NfHIW2saP+7Nq19b6iw9e+Pk+/cZ8SuusHWg8tTq+/+7dKA6v3u7LtdGfCPr+r4hjFVZataleWNnJqQ926rkeQWT07ItjtS0yg9vNLd0z2IX3m1rSc/7Ph7NrF9P/lhUNuYqjfeQyuCHd7l0MS2vbde+sMqd/bdrgv4jqbmJcm2bVVveT+izPKmKWfYxLjfOxRoUk35RxFlab1y1Gfg6Li3HY23tob0TjtTW0h+wZCtezuYmpfMbd/+kHTPkoCCIULeVB1NzUvmtm0pPFW/apf7+m5XBXy9v+OpeUmq+6RUgQM1EWXZRePk9WXE/f7VW1fLDvojynKHHyHL6rljh0zVm+tfH4e0oYOpecns9s1Uvbk6m5qXzG7bkjun6l0V8P9bH+r05jFtTvEUd98UT14PTPEcrKZR+tc6d073oOvqmmz946POv1fT2/ffPwpyWqiB/rWu46l5yfy2vbc+nGFu4pqAD4ZsvVAaRQe4ufu+yKqylW1sOzENMBYvbwqpKUgnaJLFm0Nq6HiAI8n89t0QkF7b7K5OEPFpDNhR3bLb9LYtSS+UBl11GMo1Ab9ihx3VrV+r29xTS8wXWb25rT3M7l+FeajaJuntrXSCprBtW4uivKpbKrTvRRuDSoIrZCNKb5eHtL+p8+elQtveWy+t2Ometu2agI+2A6zafMgiDY9XucMnxf3+oWBA1Vs/jCjzZmQpe/DYuLfdFYs2EvCmWL3b1vaazp8npUb73lYjfbTbPZ0g4hNtX5UKbVuSFm1wzyFWVwT8jlpbq3Z1/oOv31euptq9EWV9Bo1RWkZW3HWo3bZGIX/kQaTc4UfI8jjzEa3fZ6u0gpA3QbQ7r6nUvtmBNUNpRUgb9tF3H2zVLls7a92xA+uKgI/2phRtLqToxnMo80b07CKNQz1PJ5j0Kg7YWrItuh97KrXvJdtCqjjgjk4QXRf16D2F2rYUfaZ1N8cD3h+09WpZdI2kuq1VmAn6IttepNHzx3AO9vqWkOtOu0BsXioNKdo1N6nUvoN2eDEpkle939YbW+i72/JKWUh+FyyUdjzgN1fZUS3QkKSqthZSJGwv0B2LNA7WFFRU019wr9W7ow+xVGvfqz8h4JPZ+n22mqIcqKZa297fJFfcBdXxgC/t5MIfB2vVSCxLucOPiLsOdiio6q0fRJR5fJnKHjI+7m3Hq7TS+UaCrrFtW5ti+P5SrX1vqrRZTZ/ENtF3dyiW3353cT7go/wQGqt3q6Fie0RZn4Gj5OuVHXcdard/rGBj5Dl6ucMOl8fr/M32YvkRwV127Zfq/Z0/T0rN9l3nlz7Z72gVEIfSyuhmYFKxbUuxDV67i/MBH+WH0PaVihI0xdPGBRhyHV6k0YyV9Mkrlu8uVdv3xihDAu5D392xlA/4pqCt8ijvMNXmMZxEXQVpkzuugtSWT+qk2kbnGwpiF8vhlVRt38xQJafaRlu766J7bqq27a3VtuNXJHU04LdU2Yr272+7kSRqL9B9izQO5oZjOYhdXOtLlBrt2w2jHMQu/p1X89t20A5nnJMcPVARS3C1dZnDREzFtH0LQ59yhx3e4esCDXVadOUo+fdXSJLGnPY9TfjqbVG9Z6CxXm/cfHzL1JXHl6ljb3xBBaNntfn80kpbRwyMatNwiVgX2LmlfX/wh2tVuujeT59o6cQ7P1D2oDExvefye7+p8reeaHmvY65/Vv3Gz2nzuc0L7Xryrl+IXzK27WYfPXmj1v/rF3G/f7OxZ/xA4798c5v/tqnS1ui+CXurmDk6gq+M8kIX/vpq1e3ZHFHWu3+x0nvnxV2H/TvXK3CgNqIsZ+gEedLSO3xdWmZvlXz+spbHZS89KH9959cjtW1by+/75mfHpSxL0y57uN1wl8TtY5NQIBT9Ajs3te9Rp35HltcXfmDbKv3fb2J6v7X/uK0l3CVpykX3tRvuUnihXYDD8EmnKon77rbOm49H7oj2zwaoakjoW8XM0RG8P8ofdlXZe9Ihp9Mk7jaDXZ/iGTn/29rw7N0KNtbLX1+tspcf0pgvXtvha9Y8eaN2LvtXy+Nx59ykIUd9qcPXRHuuKdwjlu/MTe07q+9QDf3cV7X1tcclSVvf+LPGnXuzMrI7H4aUv/VXrXvq1pbHY077nobPvaDT1zUFJZ+306fBRZK5787I7qf+k+Z16T1ryj9SQ+WOlseWx9vhezp9sRtHAz6eiyTkJug4S5u3MIxy+igju69GHPctlS66T5JU+ty9GvWFK9vdg9zy2h+1/t+/bHk89Jiv6rCzru/0fQj45BNTwLusfY/54rXa+vofJdtWsLFem196SGPP/EGHr9m3/h2t/N3FLY8HzzxT479yawev+Ey0YQH3aIri1seS+9q2JE2/4vdder9P3n9B7/7q7M8KLI+mLnxIvfsXt/uaRof7bken6KOdmqve8mGrstzhExNSh4qP32lVll8yI+rXjzrlsynNhsod2vr6n9p83t51b2nVw5e3PC4Yc5SmXPq7qN4jQMAnnVhCy23tO7voMA2a9sWWx5teeEChQPuXm6z7ZJPevfNLLTf8yC+Zrunffizq4+p+2nfSMaHvjsWeNa9pyV3nRvwOJi+4R8OO/XqHr3P68JOjAZ8W5bvXbl/bqiynaFzc7x9orFflphURZb7eeTGdwpFVOExDjv5yy+MN/71LdijyW637ZFNE48jqP0JHXvt3eX0ZUb1HGtOXSSfati25s32POf37Lf/dULVT5W/9tc3n+eur9c4vz1RTzR5JUq++Q3XkdU/Jm94r6nrG8lnBHUzou6O1b/27eveXZynYdKClbOI3fqXiEy7q9LVOt21H394X5bvv37k+ssDyqFffoXG//751b8oORq6EKhw3W5YntkQd88VrpU9HK/t3btCO5f9u+beWDvDTWyX6snJ11PeeUUZOv6i3n07AJ52MGL4zN7bvglEzVDh+bsvj0udaL7YLBQNacvdXVLt9nSQpLbOPjvre08rMi+2UD9p38ol2zYQb23YsKjet1Ds/P02Bhs8uuTj+y7dq1BeujOr1TrdtRwM+mj8+6G9stVLSl5WbkHv9lr/1ZKuy/hOPj3k7OUPHa+CUk1seb/j3rySFr5O89Ndfb+kALW+aZlz9Z+UMiW0P1ulGgthF2wG6uX2POf26lv+u3vKBdn/4csS/v//YVdqz+hVJ4cVGM678Y5emX2nfyceUvrsj1VtX6+3bT5W/vrqlbOyZP9TYM77fwasipXudPf3T0YDPzez8jz94z6mZHYpyhUcH/AdqtWPpMxFlHl+Gio7seEV7e8ac/r2W/64sXaY9Hy3WB3/4rnZ/8GJL+aQL7tKALqzezMngHOFk4/NIvaJYwurm9j1g0ryIFcIb//frlv/e8N+7tPnlR1oeTzz/Fxo49WTFqlea89OYiF00fZKb23Znand8rLd+drKa9u9rKRt18tUaf+5PYtpOTnRHYbuNoz+tkfmdN5L03nmSFVnNwIHaNhtPLDY++/8UbIy81mLRrLOVkVPYpe31HXOU+h72uZbHy+/9pja98EDL45L5V2jkvEu7tO1oPie4i2VZKjagfY857bMd10/ef0E129dqx/J/a/VfbmgpL553qUrmX9Gleo7Mt7jITRIaWZD8bbs9dZ9s0ps//YIaqz9pKSs+8RJNPD/2i+M43Xc7GvAj8ix5Ovn7LY9XmbkDWpXvXfdml9/3wL5t2vDfu1uVjzzp8jaeHb0xp302pdlQtbPlvwdMnt+lxtGsJIofE9wnmu/N7e178Mwz1WfQ6PAD29YHv79Gy+/9pmSHF5L2nzRPR3yz9XtFK5qggPuURBFcbm/bbanfV643fzo/4u53w+Z8Q0dcGNsFnyTJY4UzzkmOBnxGmqUhOZ1/AIUTWl8Jq+zF6E4xO1Sw6YDevevcVnuAw2afr4JR8Z1iMXDKF5QzLPIYZM7QwzXjqj91efFHYRZT9Mkq2r13N7dvy+PR6FO/2/J4z+pXW7adM3SCZl79l7gWNkUTFHCf3ExLhVmdP8/NbftQDVW79OZt81W/Z0tLWdFRX9LUS37XpVmmoTmWMtJSOOCl6EY5RbPOalW2a+X/tGXxH2J6r8aaPXr752eo6pDTK9KzC3X41++IaVttOVCxQ/79lRFlR3zz7rjuezyKEU7SGlUQ3c/L7e172OzzlJk/OKIsI3eAjvre0/Jl5XR5u1L0nxHcx6S+u7Fmr9786RdUt2tjS9mg6V/U9G//vsuLAt0w8+r4ryuaPfhB009r8wIGKx9cqI/+epP8h6zUPFSw6YA2vfCAXvnBTO1dszji3zy+DM248vG4j98EGur07q/O0oGKbRHl5W+3ff5wtEbmO/4VoYsGZUe30M7t7duTlq4Rxy+IKJv13b8qq9/wLm9TCn82A7u+7wuHlUTRN7m9bUtSU12V3rr9FNVuW9NS1n/SPM28+i/yeLt+sVc3BLxl27ajF8v9eG9IP3yp85WV1Vs+0Os3H9/qtAtJSuuVrQGT5yt/5DRl5A6QNz1TTfsr1Vi9WxUbl6ji47cjTnVoZnnTNPPqv2jwjNPj+hvsUEhL7j5XO5f/p9W/eXyZOuk3H8d8bnCzG+ekacogQj5Z/ehlv9bs6fwn5ub2LUkfPH6dSp+7R1J49H7yA1vj3uaEfpZuPcEX93bgjJU7Q7rtteTuu/0HavXWz05W5calLWWF4+fq6P97JqaLNbXl9hPTNLbQ2b7b8YBvDNi68Bm/DkRx9sTetW/onV+cEfcqTEnK6jdc0xY+rMLxs+Pe1od//H7EKUTjz71ZG5+7p+XiNrHcSvZgPo/0yBk+9Ul3fk8QXfPEhwH9/aPorlfp1vYtSW/c+nntXfOapPDo5pgf/jfubZ47waOvTHT0dhiIw/4mWwue8Ud1WWY3tu1AY73evuM07Vv3RktZwegjdcz1zyots098dfRJj57hS+3z4KXwQru5xdFVo3DcsTr+jmURp6PFKi2zj0rmf1vH/3xFQjq/spcfjgj34hMv1tgzf6CSg1Z1Rnsr2UMdM8xDuCe5E0Z6Fe036Mb23ezg+25HezOmjlgKfzZIXn3SLR09LDn77qC/UUvuPCci3POKp+roH/w77nCXpLkjPI6Hu+SCEbwklVfbuvq5KG+erfA91fd8+LI2vfSgdn/wUqtVlYfyZvRW/shpGjzzdA2b/Y24FwY12/3BS3r7F6fLDoanH/pPOlFHff9f8njT1LS/Qs9fObplj3XCV38acRpdNO6Yl6YxfR3fB0Ocbn/Dr2Xbo/+ZuaV9N6vbs1kvXDW25fGMq/6sIUedE9c2ZxZZ+sGxTM8nu2gPsTZzQ9sOBfxacve52rXyfy1lOUMP1+dufCGq2yJH4zdf8GlIrvMB74r5saG5lg7vb2n17ug6Qcuy1H/Sieo/6UTZoaCqt3yo2h3r5K+rlL+uWh5funy98+XLylPvASOVO+zwhF+juGbbWi399ddawj27aFzEooz0PgUafvyFKv1f+PzJjc/do5IvXBn1DWZG5lsa7YJFGojf/FFeLdsefSfohvZ9sOrN70c8zhtxRNzbnD+K0bsJxvQNX9CprDI5+m47FNSye78REe59Bo/RMdf/L2Hhfnh/yxXhLrkk4KXwD3717tgvY2h5vMorntwtdxFqT2P1br3zi9NbFn9k5PTXUf/3jHxZuRHPG33Kd1T26a02G6t2aevrf4zqDkSSNH+Uhyt8GeKIgZYG9pF2deHwoxPt+1AHB3xar2z1Hjgqru0N6iNNGkjbNoFlWZo/yqP7l8V+z9+ebtt2KKQVD1ysHUv+2VKW1X+EPnfDImXmtb4gT1fNH+2enVfXzP/OHGIpL9PpWnQu2NSgd+48p+ViCB5fpmZd+3f17jei1XN7FRRp6Oe+1vJ4w3/vbnUr2bZk+aRjh7vmq0GcPJalk5J4xFq1eVXLf+cMmxj3judJo7zysPNqjGOHe5Tl8qMttm1r1SPfVvkbf24p61UwRJ/70fPqVVCUsPfJzwwffnIL16RImsfS50tcU5022batFfcvUOWGJeECy9K0hQ+p75gj233N6C9+t+V6zHW7NmrH0qc7fZ/jij2OXwEJiXV8sSfqO8y5TfWWD1r+O294fNPzPm+4fcMcmWmW67/TDx+/VptfebTlcUbeQB3zo+faHJjFY16JR2mdXX+9B7nqWzlljNfVo/g1f/uxtr/7j5bH4865SUOOPrfD12QPHqvBM05rebz+01vJtqe3TzpzXJImAdqVnWHpzMNc9XOLSmPtPh3YV97yODfOFfRnHuZRNpdeNs6Z47zq7dJR/EdP/Eili+5reZzep68+d/1zyh40JqHvk5cZzjA3cc0xeCncCS6cnqY73oz/loKJtuW1P2r9Mz9veTz02K/rsLOuj+q1Y077XsvtDavKVmr3hy+r/8QT2nzugqleFfSiAzTR2eO9Wrrd1uYqx09cidrBp8dJ8Z0iNyLP0tnj3dUBIjEKelm6cKpX9yyJ/Vh8d9r13nNa/+9fRpT5+uTrwz9Ff0/3QxUdeY5GHPfNVuULZ6S5bufVVQEvSTOHeDR7uEevb4nu4iA9Ye/aN7Tq4c/Oa+879hhNveSBDl4RKb9kuvpNmKs9Hy2WJK3/9y/bDPjpgy3NGZF8ozxEx+e1dOUsr77/QkDBJMn46oOOv1ten3KGju/SdryWdOUsr3wuODcY3WPuCI/eKQ9p+Q73NO7K0uWtyup2bYy45nysDl5X1WzOcI9mFrmv73ZfjRQexbplqn7/zg16965zFQo0SZJ6DxipWdf+XZ609Ji2M+b0z/YY96x+VZWlkTdN6O0L7wGyct5sxfkenT3elT+7NlVv/uz4e86QcTG3+2Znj/eomPsqGM2yLC2ckeaqqfqa8tUJ32beyCmRjzOlC6e6c2bKFRe6acvSbSFXTtV3l6tmeTW32J2NBInlD9r6vxcDSTVVH48ReZZ+Pi+N0XuKeLUs6Lqp+u70g2PTXDl6l1w6gpc+m6pPBUzNp5bmqfpUyDum5lPP3BEeTR+cGt+3W6fmm7m3ZpIunuZVcRS3k01mRTnSFbOYmk81xfkeXTLd/BmbS6Z7mZpPMZZl6YpZaSpK7BWTXWdkvqWLprn7N+zqX17vdEs3zUlTkaH3jO6XJf14rk85Llt5iZ4xr8Srbxzh7g4iHt+Y7NW8EnP/PrQvJ8PSj+f61C/L6Zp0j6Kc8K28e7v8ZmCuDnhJys00s6HkZUo/Oc6nwix3NxB0rzPGeXVOEi26i9Y54z064zDCPZUVZln68XE+1yyYTpR+WdKP5/iUm+n+vtu1i+wOtbfO1s2L/dpe63RN4tcvKxzug7Ld30DQ/Wzb1t8/CunJ1WYsTPrqxPBOC4edIEk7am3d/Kpfe+qdrkn8irLDs66FvZOjbSdNwEtSVYOtWxcHVJbEq4+Lcj5tIIzccYhn1wf1yMrkDvmLpnp1ssuu5gXn7a3/dIBW43RNuq44P3zIOBlG7s2SKuAlqa7J1gPLg3prq3suhBOt6YPDi0845o72vFMe0gPLAqptcromsclOly6bkaYjh5p3uAGJUdNo694lAVddCCdaxwzzaOF0r+uPuR8q6QK+2dvlIT24PKCaRqdr0rk+6eGL98wezrQlOlfVYOt3ywNasi05fppHDrF0yfQ05SXRyAbOsG1br20O6ZGVQdX5na5N53IypEunp+moJN1xTdqAl6TqBlsPrQjq7XL3juanDw5f3YnryyMWtm3rra0hPbQi6NrRfHa6dPF0r44Zyo4rYlNxwNYDy9w9mj9mmEcXT/Mm9YxrUgd8MzeO5nv7wqP2OSPo/NB1bh3NM2pHvNw6mk/2UfvBjAh4KTya/9e6oF7eFHJ0xJPlC9/v+sxx3BUOiWHbtpZut/Wfj4Nas8fZn+v4fpa+ONarmUUWO65IiIoDtp5eG9SrZSHVOxj02enSCSM9Ov0wb1ItpOuIMQHfrClo6+2tIS3aGNL6fT33pxXnW/rCKI+OHe5RRpoZjQPus6UqpOc3hrR4c0gNPXSrhsy08OVHTxrl0fC85B/VwJ0aArbe3BLScxtDKqvsub57TF9L80d5dPQwj9INu6SycQF/sNKKcGf4+paQmrrh7COfJ3ycZv5oj0YXMKJBzzngD09vPrcxpPLq7vkJD8sNd3xzRnjUy0fbRs+wbVsbKmwt2hDSW1tD8nfDEqt0rzR7eHintaTA3J1WowO+Wb3f1oZ9tkorbZVWhLSpwtYndbFvpzBLGlVgaWS+RyUFlkb3tdQnyU6bgFls21Z5ja3Sis/+V1Zlx7xDm+6VivMslRR89r+hOey0wln7mz7tuytslVaGVFpha28XLpgzoLdUUuDRyAJLJfnhvjsrBXZaUyLg21LbaGtTZTj0qxtsLX79La1eu14eX4ZC/kaNLinWSSfOVU6GpZH54Q4vmVdTInUEQ7a21YTb97YaW1W1B/TY43+SJ72XJCnUdEDf+sZ5ysvupSE54fY9JMeS10P7hvtVN4Tb9qZKWzWNtp5/abE2lJa19N2HjxujubOPUW5mOMxH5lvKTtG+O2UD/lDf/e53dffdd7c8vvDCC/XII484WCMgMfbs2aP+/ftHlO3evVv9+vVzqEZA4ixYsECPPvpoy+NrrrlGd911l4M1cg9zDz4AAJDCCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABlpWV/cAACAASURBVCLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAA6U5XYGe5vf72ywPBoMRj0OhULvP9fl8Ca8XEK9QKNSqHUttt3m/399mudfrlcfDfj/cp73+OBQKRTwOBoP03Z+ybNu2na5ET/jggw/01a9+VWvWrIl7W6NHj9Yf//hHzZo1KwE1A+ITDAZ1zTXX6JFHHlF9fX1c28rKytKCBQt09913y+v1JqiGQNctWbJE559/vjZs2BD3tsaPH68nnnhCkyZNSkDN3C9lAv64447T4sWLE7a9qVOnasWKFQnbHtBVL7zwgk466aSEb3PevHkJ3SbQFdOmTdPKlSsTtr3jjjtOr7zySsK252YpMxfX0NDg6u0BXdUdbZH2Dbc4cOBAQreXSm07ZQL+vPPOc/X2gK468cQT1b9//4Rtb8CAATrhhBMStj0gHvTdXZcyU/QNDQ0aNWqUtm/fHve2CgoKtHnzZmVnZyegZkD87rrrLl177bUJ29Y111yTkG0B8aqpqVFxcbEqKiri3lZRUZFKS0uVkZGRgJq5X8qM4DMzM/XDH/4wIdu67rrrCHe4ysKFCxMyih8wYIAuvfTSBNQISIycnJyE7bxef/31KRPuUgqN4KXEjOIZvcOtEjGKZ/QON0rEKD7VRu9SCo3gpcSM4hm9w63iHcUzeodbJWIUn2qjdynFRvBSfKN4Ru9wu3hG8Yze4WbxjOJTcfQupdgIXopvFM/oHW7X1VE8o3e4XTyj+FQcvUspOIKXujaKZ/SOZNGVUTyjdySDroziU3X0LqXgCF7q2iie0TuSRayjeEbvSBZdGcWn6uhdStERvBTbKL5v374qKysj4JE0YhnFM3pHMqmtrdWIESOiGsWn8uhdStERvBTbKJ7RO5JNtKN4Ru9INtnZ2bruuuuiem4qj96lFB7BS9GN4hm9I1lFM4pn9I5kFM0oPtVH71IKj+Cl6EbxjN6RrDobxTN6R7KKZhSf6qN3KcVH8FLHo3hG70h2HY3iGb0jmXU0imf0HpbSI3ip41E8o3cku4ULF6qwsLBVeb9+/Ri9I6l1NIpn9B6W8iN4KTyK79evn/bv399SlpGRoT179hDwSHq33nqrbrrppoiy2267TTfccINDNQISo7a2Vv369VNjY2NLWZ8+fbR3714CXozgJYVH8RdccEFE2VlnnUW4wwjXXnut0tLSWh6npaUxNQ8jZGdn66yzzooou+CCCwj3TzGC/1QoFNLZZ5+t119/XTNmzNB///vfiE4RSGbPPfecLr74YknSQw89pC984QsO1whIjEAgoFNOOUXLly/X7Nmz9dRTT8njYewqEfAAABgpZYeotY22Sittbaq0VXXAlj8kNQWkQEhK80g+r5TulXIyLI0ssFSSbyk303K62kCngiFb22pslVaE/78xKDUFJX8w/O/NbTvDKw3JsVRSYGlIjiWvh/YN96tu+LTvrrBV02i3tO3mvjs9TfJ5pLxelkbmh/vu7IzUbNspEfD1flvr94UbRGllSKUVtnbXxb6dwiyppMBSSb5HIwssjelrqU96ajYcuINt29paHQ7z5k6vrCrc6cUi3SsV5322M1tSYGlYriXLon3DOfubWvfde+tj307/3q377iyf+W3b6Cn60oqQFm0M6Y0toZg7vGj4PNLRwzyaP8qjMX3pDNFz6v22Xtscbt/l1d3zEx6aa2n+KI/mjPCkRGcId7DtcKgv2hjSW1tDCoQS/x7pXunY4eG+u6TA3OP1xgV8U9DWW1vDHd+GfT33pxXnhzvDY4d7lJlGZ4jusaUq3LZf2xxSQ6Bn3jMzTZozItwZDs8ztzOEsxoCtt7YEm7fZZU913eP7hvuu48Z5lG616y+25iAr26w9a91Qb28KaTaJufqkeWTjiv26MxxXhX0MquxwBm2bWvpdlv/+TioNXuc/bmO72fpi2O9mlnEjBUSo+KArafXBvVqWUj1fufqkZ0unTDSo9MP8xqz3sqIgH+7PKQHlwdU09j5c3tKb5+0YKpXc0Z46AjRZVUNtn63PKAl29z1Mz1yiKVLpqcpz5COED3Ptm0t3hzSoyuDqnMw2A+VkyFdOj1NRw1N/tmqpA746gZbD60I6u3ybjhIkyDTB1taOCON0TxiYtu23twa0kMrgtrv4IxUR7LTpYuneXXMMHZiEZuKA7YeWBbQ8h3ujZ+jh3p08bTkHs0nbcC7cdTeHkbziIVbR+3tYTSPaLl11N6eZB/NJ13A1zXZun+Zu0ft7Zk+2NIVs9KUk6LnZKJz75SHdP+ygGtH7e3JTpcWzkjejhDdr6bR1r1L3D1qb8/RQz26bIZXvZPstOikCviqBlu3Lg6orCppqtxKUY7047k+FWYlV0NB93t2fVCPrOyG8zl70IKpXp0yxut0NeAye+tt3bzYr+01Ttek64rzLd04J7lmqpIm4PfUhRvIjlqnaxK/flnST47zaVB28jQUdB/btvX3j0J6cnVyh3uzrxzu1ZcmcDgKYTtqbd38ql97unCBGrcpyv50gNY7Odp2UgT83jpbN7xsRgNplpcp3XaCT4MJ+ZT3lw8C+sea5Dvk1JFzxnv0tUkpcaFMdGBHra0fvexXVYPTNUmcflnST09IjpB3/QGzqobwyN2kcJekqgbp5lf92lvv+v0rdKNn1gaNC3dJ+seakJ5ZZ8aMBLpmb72tn7xqVrhL0p566ebFflU3uL/vdnXA1zXZuvW1gLYbMC3fluaGUtPo/oaCxHuxNKjH3zc3BB9fFdSLpeb+fWhfTWN4YNaV68Yng+210i2vBVTX5O6+29UB/9CKYI9estAJ22uke5cElARHSpBAmypDenC5+eH34PKgyirNm6FA+2w7vFo+mRfURaOs0tbDK9z9G3btQbKl20J6fUt8HcPav9+idf/8aURZ/uhZmvOTV2V5ur7S17ZtvX37Kdr94csR5SXzv61JF9wV8/aW7wjfOGRuMauPU4E/aOueJUEFu7hPZ4eCeu2mOaosXRZRftjZN2rcOT9KQA2lxpo9evl7U9VYszui/Mjv/VODpp4S9XaCtnTPkqB+Ps+Sz7DrfKNtizeH4joVrm7PZr1w1dgE1qhzh539I40758aYX/falpCOGhbSzCJ3jpVdWavaRlsPLI//Thpjz/yhcocfEVFWuWGJNvz37ri2W/bi71qFe59BozXhqz9t5xWde2RlUBUHGMWngqfWBLUljlM9LY9X0y5/RN70XhHlHz9zh6rK3ou3epKkVY9c2Srch8/9Zkzh3mxzla2n1rh7pIPEqDhg69EkP9UzVg8sC6jWpYdZXRnwj6wMJmRhhifNp2mXPSJPWnpE+dp/3KKa8jVd2ub+T0q1+i/XR5RZHq+mXda6w41FnT/cUJiqN9umypCeSsCiuuzBYzX+K7dGlNlBv5b/9kIF/fFd3rH8zSe0Y+nTEWVZhcM18fxfdnmbT60JMVVvONsOX342Ga5Ql0hVDXLtTo3rpugTMTV/sNzhE3XY2T/Smr/e1FIW8jdqxf0LNOfWN+TxRv8R2KGQVtx/kYKNdRHlo069RgWjZ8VdV6bqzRbv1PyhSuZfoZ3L/6O9a15rKavdtkZr/vZjTfz6HV3a5oGKHXr/99dEFlqWpi58UL6snC7Xlal688U7Nd/M401X7vBJCahRa3WfbFKgYX+r8qzCYXFt161T9a46D7620dbVzyX+tIr2jlmOO+cmHXb2DVFvZ/1/7tRHh4zec4ZO0HE/e7fVLEFX9fZJvz7Zx81pDPTkhwH97aPEjmLr92zRy/83TYEDB51qYnl07E0vqfCwY2Le3ts/P02frHo+oqxk/hWadMGd8VZVknTuBI++MtF14wrEqeKArav/53f16H3fx2/rjVs/LzsYWcmio76kmVf9Ke7t52VKv/6CT9kuuhS5q3Y3nl2fmKn5Q1ker6Ze9pA8vsyI8nVP366qze9HtY2a8jVa+/ebI7frbfsQQDzq/NLTa9053YOuq2209fS6xE9RZ/Ubronn/yqy0A5p5f0L2hypdKTs5UdahXufwWM04au3xVvNFk+vC7n2eCW67um17r55zIGK7Vpy91dahXvOsImaeumDCXmPqoZwhrmJawI+ELL1Qmn3HaPLKRqn8ef+JKLMDvq14rcXKhTo+M4eoWBAK+6/UKFDjm2OPeMHyiuekuiq6tWykBoCdIImeaUsJH83/fZHHPdNDZx6ckRZ3e4yffin/4t6G3V7Nmv1nyOfH15b8mhca0sO5Q+G2zfM0RCwXf2dBv2NWnLXl9VY/UlEua9PgY689u9Ky8hK2Hu9WBpSIOSevts1Ab90m93tVzwadfLVKhh7dERZTflqrf3Hre28Iuzjp29vtTo5r3iKxp75g4TXUZLq/dKbCVyHAGeFbFvPb+zePfspF9+v9D59I8o2v/xwqxF5W2zb1soHLomc5pc05vTvq2DUjITWU5Ke3xhUyD1HBhGnN7aEVO/i0fuqR65sdXjW8ng186o/qXf/4oS+V2WDtHS7e9q2awJ+UTd3gJJkeTyatvAheQ/ZY9vwnztVsXFZm6+p3LRSHz8TuWDJ48vQtMsejWmBXqye2xhiRb0h3t9la1dss+Uxy8wbqMkL7mlVvvLBhWraX9nha0sX3RuxUE+ScocfocPOin59Six27pc+2EXbNoFt21q0wb2DkdJF92nra39oVT7hq7ep/8QTuuU9F21wzzS9KwK+vNrW6t0984PvM3BUq/PV7VBQK+5foGBT5BRC0N+oFb+9UHYw8pz8cefcpJyh47u1nmWVtjZU0AmaoCd2XiWp6MizNeToL0eUNVTu0Pu//067r6nduV5rnoy8wIfHl6Hp335MnjRft9RT6rnPBN1r/T7btbfv3rvmdX34p++3Kh9y9Jc1+tTvdtv7rt5ta1u1Oz4TVwR8T//YR37+MhWOnxtRtn/HxxGn0knSmr/9WLXb10aUFYw+UqNPPeQ0om7i5j1jRGd3na3lPThld8S3fq3M/MERZdveelLbl/yz1XPDO7YXKdh0IKI8vAM7oVvruWy7rT117ugE0XWLNrqzj6rfV64lv/5aq8FZ7vAjNOWSB7r9/d2yA+t4wDcGbC3u4QUalmVp2sIHldYrO6J843P3aO+6tyRJe9e9pY3P/jri370ZWZp22cNxXeY2Fm9tDWm/y29mgI69vCmonvwG0/vka8ol97cqX/XIlWqoilxktP4/d6lyw5KIsoKxR3fr6KaZrfBng+S1v8nWW1vdF/DBpgYtufNcNdXsiShP79NXs679W0IX1bVn8eaQGl2wUNrxgN9cZetA/FeljVlWv+E6/Os/jyy0Q1r5wEVqrN6tlQ9cJNmRjXfCV25Tn0Gje6yO/pC0YZ/zjcRkoVBI1113nYYNG6ZTTz1Vb7zxRkK3/1EPHXo62MDJ8zXihAURZU21e/XeQ5e1PK4p/0jr/nFLxHO8Gb0/3YHtmW7Bic8m1SxZskTTpk3T+PHj9bOf/Uy1tYm7Nef6fbYC7st3vffQZaoqWxlRZnm8mnH1n9S734geqUO9P5xtTnM84EsdPM5cfMIC9Z80L6Ks7pNNevkHM1T3yaaI8sLxczXypMt7snqSnP18UsEzzzyjO++8U+Xl5Xr22Wc1e/ZsnXDCCQkJ+pBtO3Y3xInn/UJZ/UdElO1a+ay2LP6DQoHwJW0PPT308K/frj4DSnqsjpsqbVbTd7Pzzz9fK1eu1Nq1a3XDDTdoxIgRuv322xMS9Jtc2DdtePb/qfzNv7Qqn/C1n6n/4cf3aF02ueBOqM4HvMMfwtRLfydf77yIssaqXRGP03pla9rCB2VZPX+Fok1cv7tbrV69ulXZK6+8kpCg31krR2anJCkts4+mLXxYsiJ/4h88fp3ef+wqVW9eFVHef9KJGjnv0p6sog4EpF2JG1DiEH6/Xxs2bIgoq6io0PXXX5+QoC91Wd+0e/Urra40KklDj/mqRp/S/kLT7uKGwZnzAe/wh9CroEiTvtHxZTgP//rPldVveA/VKNJGFzSSVBVv0G+scLYDLBx3rEadfFVEWeBAjTa/8mhEma93nqZe8ruerFoLpz+jVJWIoHe67z5Y3Z7NWvbr82SHItd15I6Y3OaalJ7ghs/H0YBvDNjaVuP8hzBs9nnqP+nENv9twOSTVHzI8cyetLdequHSno7qatC7YYpu/Lk3K7toXIfPmXTB3erVd0gP1SiS0zN4qa6rQV/dYGtvfTdXLkqBxnotufNLatq/L6I8PbtQR373bwm9EmMsymtsxxfaORrwm6tsueGqfo01e1W95cM2/23wzLN6uDatuWFPELEHvRu+N296pqZd/qisdi7KNHjmGRp27Nd6uFafceNx3FQUa9C7acds5e8uUfWWDyLKLG+aZl79Z8dmXiUpZDu/0M7RgHfDCEeS3n/0qlbXKW720ZM/UuMhp1v0NLd8TgiLJuhtBxfYHSp/5FSNPaP1ZZXTsws1ecG9DtToM5sqba7Y6CLRBr1bdszW//tX2v7O31uVH/6129Vvwtyer9AhnO67HQ346gbnG8m2t/+m7Uueavffm2r26P1Hr+7BGrXGFL07dRT0/pBzC+zaUjh+dquynCHjlZHTz4HafOZAQK481SrVdRb0buiTPvngxVYXJ5Okocd+vdXaE6fUNHb+nO7k6I2Zmxy+zkVD1S6teqx1ePuycuWvr255vH3JU9r+7lMqOvLsnqxei+dfWqw3733Ckfc23Xvvvdf5kzrxyiuv6JVXXtHw4cP1m9/8Rqeddlq33TnORJdefoWsQDffaSoFhULx7zk1B/3NN9+ss846S48//rjS0tIc77v3f1KqZb85v9WiurziqZpy0W8dqlVrTUFnd4QcDXi/w3vu7z10ufz7KyLKik+8WP0nzdOSu86NKF/12NUqHD/bkRHPhtIyrXz00c6fCEdt2bJFZ5xxhtauXat+w8Y4XZ2k8Ze/PtXuITK4Q2Njo5544glVVFRo0aJFju7ABhr2a8mdX5K/LvImSuk5/TTru3+VNz3ToZq15vSOkKNT9E5OzW157XHtWvlsRFlW/xE6/Ot3aPCM01vdtKOpZo9WPerMtI/Hl+HI+yJ2tm3r2WefZdo5Bp60dKergCgtWxa+66aT7XvF/ReppvyjiLKWRXWFwxyqVduc7gccDfg0h979wL5t+vDx6yILLUvTLn1IaZl9JElHfOv/KSNvYMRTdiz5p7a984+eqmaLkN/hAzmImsfj0VlnnSWf41eYSB6hAO07WcyeHV7L4VTf/fHTd2jH0qdblU887+fqN36OAzXqmFOfU8v7O/nm6T1zz5ZWVj54acQxdkkqOenyiIVI6X0KNHnBPVpy55cinvf+Y1er3/jZysjt3yN1laTDx43RnGt65g52qebdd9/VO++8k5BtTZgwQffdd59GjBihOm4SFLWFF18oT5CQT7RgMKjf/OY3CdlWnz59dMEFF7RsL92B5Nj13nNa8/ebW5UPm32+SuZf0fMVikKGQxnXzNGAd2KUU/big9r9wUsRZb0HjtL4r9zW6rmDp5+moZ/7WsS1jZtq92rVo1dp1jVPdntdm82dfYy+NcV9e6cmuOWWW+IO+DPPPFM33XSTJk+e3FLm1M5rMvrlHT+Tz9vzl4E2nd/vjzvgi4qKdP3112vBggXKyPjsUGFP9921O9dr+b0XtLoBWN7IaY6f6tkRp9u1oxMI+b169o+v212m1X/5YWSh5dG0yx5u9xaCky64S5l5gyLKdix9WtvaOPeyu+Rm0vm50Zlnnqn33ntP//znPyPCXQpPzWX5HKpYEuntc34aE60VFRXpvvvuU2lpqS6//PKIcJekvB7su/0HasOL6g6Zdc3I6e+6RXWHynO4ao7+tEbm91wjsW1bKx+4RIGG/RHlo065Wn3HHNXu69L75GvyRa33EN9/7DtqrN6d8Hq2paQHPyd0rqNgb2ZZVo+272Q1Mt9y5CZOaFtnwd6sp9q2bdta8dtvqXb7uohyy5ummd/5i7L6Du2RenSV032AowE/PM9ST81glD53j/aufT2iLLvoMI3/0k86fe2gaadq6LHnRZSFp+qvTGQV2+V0I0FYNMF+sJICvrfO8Bm5Q7TB3qynBh3rnrpNO5f/p1X5xPN/qcJxx/ZIHbrKa4UzzkmOBny619LQ3O7/AGp3rm91xSPL49W0yx6Jenpn0gV3KjN/cETZjqXPaNvbf0tYPdsyoLeUnUEn6KRYg70ZMy+dG0nAOyrWYG+WnWGpf+/urdvO5f/Run/+tFX5sDnfUMlJl3fvmyfAsFxL6al8DF7q/j14OxTUivsvUrDpQET56NOuU37J9Ki3k947T1Muuq9V+fuPfUcNVd13kY6SAse/IqN1ND3c1WBvxnfXuVH5fEZO6GqwH6w7++7a7eu0/LcXSofcpyC/ZIarF9UdzA2zU47/urp7lLP+P3epcsOSiLKcoYdr3Nk/inlbA6eerGFzvhFR1rR/X7dO1TPC6V5tBXe8wd5sYB8W2nWkt08a0MfpWpjL5/Np/PjxEWWJCPZmJd20c+avr9a7d56jwIGaiPKM3AGadc2T8ibJhb/cEPCOniYnde+HUFP+kdb945aIMsvr07TLHuny1bMmnv9L7f7wZTVUbG8p27nsXyp/668aesyXO3hl1zDN271OPfVU3XLLLfrb3/6mI444Qtddd11coX6w5oV2q3dzTnxbWGDX/Z544gl95zvfUUNDg84777xWp7vFo7sGH+///rvav3NDq3LbDumdX57ZLe8pScffsSyh23PD2inHA35EnqU+6dL+psRuNxTwa8X9CxQKRG547Bk/UF5x1zvw9N55mnLxb/XOz0+PKP/g99eo34S5yswb0OVtt3ovrzS6r/ONxGSWZenGG2/UjTfe2C3bP7y/R6t3c+eZthw+wPEJRONNmjRJr7zySrdse0xfS+nexF9v/cDerW2WN9XsUZPDt+6OVp/0cLY5zfFfmM9r6bjixFfj42fuUFVZ5J3C8oqnaOyZre+LHauBk+dr+NxvRpQ17d+nVY8k9mpKs4d7lOVzvpGg604s8cjDV9iK15JOGOl494M4ZPksHTuc77Atxxd7HL/IjeSCgJekk0Yl9rJfVWXv6eNn7ogo86Slh6fmvYmZtJh4/i/Uq2BIRNnO5f9W+VuJu8LdSaNc8fUgDgW9LM0a4vwP3W1mDfGooIcvdIXEm08f1aZEZ1pXWbZtu+IA4S2L/Vq1PKFQdQAABr5JREFUyxVVcYUxfS3dMY8VWib48JOQfvxqwOlquMotx6UxRW+I/3vRrw376LubTR5o6aa57ui7XfMLm++SPR63YM/YHIf3t1SU43Qt3GNIjjShP6N3U9BXRZo/2j1Z5ppvZtpgS4VtXw4+5WSnS0cPc81XgzhZlsUO7EHmj/Kyet4gxwzzqE/XTkoyTmGWNG2Qe9q2a1LE67H0+RI6QSm8+MjpKyAhseaO8CjT8XNWnJeZJs0Z4ZpuBwmQ7rV0IgsmJUmfL/HK66JVta76Vk4e40n5UXxOhnT6YezomKZ3uqVzJvC9fmmCV73T3dMBIjFOP8yrnOS4/ky3KcwKZ5ibuKo2WT5Ll89I7WHOpdPTuD2soU4f69FoF1zdyimj+1o6bayruhwkSG6mpUump3bfffnMNNed1uy6X9vkQR7NS9HpnmOGeXTU0NT821OB12PpillpKXn/c59HunJmmqumL5FYRw/16OgU7b/mlXg0eaD7/nb31UjSBVO8KTdVn5MhXTSVKVzTDc219JWJqfc9f2WiV0N64M6RcNbF01Jvqr4wS7pgsjt/064M+FScqmdqPnWk2lT96AKm5lNFKk7Vu3Fqvplrf3WpNFXP1HxqSaWpep9HunIWU/OpJJWm6t06Nd/MvTWTdOFUrw4rNLtjKM63tHC6O6d30H2G5lq6apZXJrduS9KVs5iaT0WXzfCq2AV3U+tOhxVaunCKu/tuVwd8RpqlG2anqdgFd+XpDkXZ0k1z0jhtKEV9brhXlxi8c3fpdK8+N9zcvw/t651u6cY5aRqc7XRNukdxXjibMtLc3Xe7OuClTxvK3DQVGdZQ+mVJP57r47h7ijtplFfnH2FeCJ5/hFef5+p9KS0v09JP5vrUz7AF00XZ0o1zk2Ng5pqbzXSmqsHWrYsDKqtKiup2qCgnHO6FWe5vIOgZ//04qEffM+O+8QumenXKGMIdYXvrbN282K/ttU7XJH7FeeEBZ16SDMySJuAlqa7J1k9fD2jd3qSpcisj88NTV4zccahXy4K6d0lQydq6PZb07ZleHVdMuCNSdYOtW14LqKwyWVt3+Jj7DbOTY+TeLKkCXpIaArYeWxnUi5tCTlclZnOGe3TRNC7Vifat2hnSvUsDqjjgdE1iU9BLumJmmiYPcv1RPzikrsnWwyuCem1L8vXd80o8unCK1/XH3A+VdAHfbNXOkH67LKC99U7XpHN5mdLCGWmaWUTnh87VNdl69L2gXi1Ljo7w+GKPvjWFHVdEZ+n2kB5YFlBVg9M16VxhVvg8dzefCteRpA14Sar32/rDe+4ezc8Z7tGFU73KzqDzQ2xW7Ajp/mXuHc0X9JIum5GmaYOTs/ODc2obbT2yMqjXXTyan1fi0QWTva69iE00kjrgm7lxNM+oHYng1tE8o3YkwtJtIT2w3F2j+WQftR/MiICXwqP5/60P6YXSoKNBn58Z3vM7ZQyjdiTOh5+E9J+Pg1qxw3ZsEZ4ladpgS18c69XEAcnf+cEdahttPbs+qBdLQ6p0MOgLs8L3cz95jCepR+0HMybgmwVDtlbstLVoQ1CrdvXcn3Z4f0vzR3s1s8hSGpflRDfZvd/WC6VBvbQppJrGnnnPnAzpxJEefb7Eq/59aNvoHoGQraXbw3336t0913dPGWRp/iivpg6yjLuksnEBf7Cdtbae3xjUK2Uh7W9K/PazfNLcER7NH8XlONGz/EFb75SHtGhjqNtOGz2s0NL8UeH7JPi8tG/0nPLqcN+9eHNI9f7Eb79Pevgw00mjvBqUbW7bNjrgm/mDtjZX2dpUaWtjha1NFba2VtsKxvCXeyxpaI6lkoLw/0bmWyrOt5ROxweH7a23VVphq7QipNLK8H/HOrrPyVC4bedbKinwqKTA4kJMcFxjoHXfXV5jKxRD3+21pGG5kX33iDwrJXZaUyLg29IUtLXl04ZT1RDeCWgMSoGQlOaR0r1SutdSToZaGkSynQOJ1GTbtvbVS6WVtrbV2GoM2GoKSk2fXigv3LbD93oYkhMO9b5ZkmXRvuF+B4d+TWO4L286qO/O8Eo+r6W8zHDfPTwvdQdiKRvwAACYjKWwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAMRMADAGAgAh4AAAMR8AAAGIiABwDAQAQ8AAAGIuABADAQAQ8AgIEIeAAADETAAwBgIAIeAAADEfAAABiIgAcAwEAEPAAABiLgAQAwEAEPAICBCHgAAAxEwAMAYCACHgAAAxHwAAAYiIAHAMBABDwAAAYi4AEAMBABDwCAgQh4AAAM9P8BWpOoq2f68WwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 660x660 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(2.2, 2.2))\n",
    "G2.plot(ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And save it as a dot file so that we can reuse it later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "G2.save_model('dags/Primer.SCM.1.5.3.dot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this is a stochastic model, we have to use random numbers. We'll assume that the exogenous variables are normally distributed and generate 100,000 values for each of them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100_000\n",
    "Ux = np.random.normal(size=N) \n",
    "Uy = np.random.normal(size=N)\n",
    "Uz = np.random.normal(size=N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And define the 3 functions that determine the values of X, Y and Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "fX = lambda ux: ux\n",
    "fY = lambda x, uy: 1/3*x+uy\n",
    "fZ = lambda y, uz: 1/16*y+uz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these two ingredients at hand, we can easily build up the DataFrame that describes our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = fX(Ux)\n",
    "Y = fY(X, Uy)\n",
    "Z = fZ(Y, Uz)\n",
    "\n",
    "model = pd.DataFrame({'X': X, 'Y': Y, 'Z': Z})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And inspect it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>Z</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.536569</td>\n",
       "      <td>-1.402917</td>\n",
       "      <td>1.009247</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.468823</td>\n",
       "      <td>-0.351639</td>\n",
       "      <td>-0.021071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.180558</td>\n",
       "      <td>-0.749858</td>\n",
       "      <td>-1.535138</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.562355</td>\n",
       "      <td>-2.155580</td>\n",
       "      <td>1.335706</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.101650</td>\n",
       "      <td>-0.735451</td>\n",
       "      <td>0.032497</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99995</th>\n",
       "      <td>0.978800</td>\n",
       "      <td>-0.137405</td>\n",
       "      <td>0.742942</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99996</th>\n",
       "      <td>1.170064</td>\n",
       "      <td>-1.348096</td>\n",
       "      <td>1.796047</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99997</th>\n",
       "      <td>0.763898</td>\n",
       "      <td>0.180325</td>\n",
       "      <td>-1.051737</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99998</th>\n",
       "      <td>-0.491289</td>\n",
       "      <td>0.259639</td>\n",
       "      <td>1.337458</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99999</th>\n",
       "      <td>-0.762940</td>\n",
       "      <td>-1.345078</td>\n",
       "      <td>0.302126</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>100000 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "              X         Y         Z\n",
       "0     -0.536569 -1.402917  1.009247\n",
       "1      2.468823 -0.351639 -0.021071\n",
       "2      0.180558 -0.749858 -1.535138\n",
       "3     -0.562355 -2.155580  1.335706\n",
       "4     -1.101650 -0.735451  0.032497\n",
       "...         ...       ...       ...\n",
       "99995  0.978800 -0.137405  0.742942\n",
       "99996  1.170064 -1.348096  1.796047\n",
       "99997  0.763898  0.180325 -1.051737\n",
       "99998 -0.491289  0.259639  1.337458\n",
       "99999 -0.762940 -1.345078  0.302126\n",
       "\n",
       "[100000 rows x 3 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we observe that Y=3, then the expected (mean) value of Z is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.18488105416791895\n"
     ]
    }
   ],
   "source": [
    "Z_Y3 = np.mean(fZ(3, Uz))\n",
    "print(Z_Y3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is close to the analytically expected value of 3/16=0.1875"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way of calculating this value would be through a direct fit using the numerical values in our DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>Z</td>        <th>  R-squared:         </th>  <td>   0.005</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   0.005</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>   489.2</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sat, 25 Jul 2020</td> <th>  Prob (F-statistic):</th>  <td>3.96e-108</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>17:37:39</td>     <th>  Log-Likelihood:    </th> <td>-1.4195e+05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>100000</td>      <th>  AIC:               </th>  <td>2.839e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td> 99998</td>      <th>  BIC:               </th>  <td>2.839e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -0.0026</td> <td>    0.003</td> <td>   -0.823</td> <td> 0.410</td> <td>   -0.009</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Y</th>         <td>    0.0666</td> <td>    0.003</td> <td>   22.117</td> <td> 0.000</td> <td>    0.061</td> <td>    0.072</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.435</td> <th>  Durbin-Watson:     </th> <td>   2.008</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.488</td> <th>  Jarque-Bera (JB):  </th> <td>   1.439</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.003</td> <th>  Prob(JB):          </th> <td>   0.487</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.983</td> <th>  Cond. No.          </th> <td>    1.05</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      Z   R-squared:                       0.005\n",
       "Model:                            OLS   Adj. R-squared:                  0.005\n",
       "Method:                 Least Squares   F-statistic:                     489.2\n",
       "Date:                Sat, 25 Jul 2020   Prob (F-statistic):          3.96e-108\n",
       "Time:                        17:37:39   Log-Likelihood:            -1.4195e+05\n",
       "No. Observations:              100000   AIC:                         2.839e+05\n",
       "Df Residuals:                   99998   BIC:                         2.839e+05\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.0026      0.003     -0.823      0.410      -0.009       0.004\n",
       "Y              0.0666      0.003     22.117      0.000       0.061       0.072\n",
       "==============================================================================\n",
       "Omnibus:                        1.435   Durbin-Watson:                   2.008\n",
       "Prob(Omnibus):                  0.488   Jarque-Bera (JB):                1.439\n",
       "Skew:                          -0.003   Prob(JB):                        0.487\n",
       "Kurtosis:                       2.983   Cond. No.                         1.05\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = smf.ols('Z ~ Y', model).fit()\n",
    "ols.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.197132\n",
       "dtype: float64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols.predict({'Y': 3})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the added advantage that this approach takes into account also the values of the exogenous (error) variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If, on the other hand, we observe that $X=3$ then our calculation requires 2 steps, as this information must *flow* from X to Y and finally to Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.05970934593248185\n"
     ]
    }
   ],
   "source": [
    "Y_X3 = fY(3, Uy)\n",
    "Z_X3 = np.mean(fZ(Y_X3, Uz))\n",
    "print(Z_X3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which again is close to the expected analytical value of 1/16=0.0625"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or a simple fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>Z</td>        <th>  R-squared:         </th>  <td>   0.001</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   0.001</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>   51.92</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sat, 25 Jul 2020</td> <th>  Prob (F-statistic):</th>  <td>5.83e-13</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>17:37:39</td>     <th>  Log-Likelihood:    </th> <td>-1.4217e+05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>100000</td>      <th>  AIC:               </th>  <td>2.843e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td> 99998</td>      <th>  BIC:               </th>  <td>2.844e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -0.0028</td> <td>    0.003</td> <td>   -0.879</td> <td> 0.380</td> <td>   -0.009</td> <td>    0.003</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X</th>         <td>    0.0228</td> <td>    0.003</td> <td>    7.205</td> <td> 0.000</td> <td>    0.017</td> <td>    0.029</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.508</td> <th>  Durbin-Watson:     </th> <td>   2.009</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.470</td> <th>  Jarque-Bera (JB):  </th> <td>   1.511</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.003</td> <th>  Prob(JB):          </th> <td>   0.470</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.982</td> <th>  Cond. No.          </th> <td>    1.00</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      Z   R-squared:                       0.001\n",
       "Model:                            OLS   Adj. R-squared:                  0.001\n",
       "Method:                 Least Squares   F-statistic:                     51.92\n",
       "Date:                Sat, 25 Jul 2020   Prob (F-statistic):           5.83e-13\n",
       "Time:                        17:37:39   Log-Likelihood:            -1.4217e+05\n",
       "No. Observations:              100000   AIC:                         2.843e+05\n",
       "Df Residuals:                   99998   BIC:                         2.844e+05\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.0028      0.003     -0.879      0.380      -0.009       0.003\n",
       "X              0.0228      0.003      7.205      0.000       0.017       0.029\n",
       "==============================================================================\n",
       "Omnibus:                        1.508   Durbin-Watson:                   2.009\n",
       "Prob(Omnibus):                  0.470   Jarque-Bera (JB):                1.511\n",
       "Skew:                          -0.003   Prob(JB):                        0.470\n",
       "Kurtosis:                       2.982   Cond. No.                         1.00\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = smf.ols('Z ~ X', model).fit()\n",
    "ols.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.065724\n",
       "dtype: float64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols.predict({'X': 3})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also run the computation in the opposite direction (contrary to the direction of causality). If we want to know the expected value of $X$ given that $Y=2$. In this case we must fit a linear regression model to obtain a the dependency of X as a function of Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>X</td>        <th>  R-squared:         </th>  <td>   0.101</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   0.101</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>1.119e+04</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sat, 25 Jul 2020</td> <th>  Prob (F-statistic):</th>   <td>  0.00</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>17:37:39</td>     <th>  Log-Likelihood:    </th> <td>-1.3664e+05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>100000</td>      <th>  AIC:               </th>  <td>2.733e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td> 99998</td>      <th>  BIC:               </th>  <td>2.733e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -0.0011</td> <td>    0.003</td> <td>   -0.352</td> <td> 0.725</td> <td>   -0.007</td> <td>    0.005</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Y</th>         <td>    0.3019</td> <td>    0.003</td> <td>  105.766</td> <td> 0.000</td> <td>    0.296</td> <td>    0.308</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 0.988</td> <th>  Durbin-Watson:     </th> <td>   2.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.610</td> <th>  Jarque-Bera (JB):  </th> <td>   0.999</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.006</td> <th>  Prob(JB):          </th> <td>   0.607</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.991</td> <th>  Cond. No.          </th> <td>    1.05</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      X   R-squared:                       0.101\n",
       "Model:                            OLS   Adj. R-squared:                  0.101\n",
       "Method:                 Least Squares   F-statistic:                 1.119e+04\n",
       "Date:                Sat, 25 Jul 2020   Prob (F-statistic):               0.00\n",
       "Time:                        17:37:39   Log-Likelihood:            -1.3664e+05\n",
       "No. Observations:              100000   AIC:                         2.733e+05\n",
       "Df Residuals:                   99998   BIC:                         2.733e+05\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.0011      0.003     -0.352      0.725      -0.007       0.005\n",
       "Y              0.3019      0.003    105.766      0.000       0.296       0.308\n",
       "==============================================================================\n",
       "Omnibus:                        0.988   Durbin-Watson:                   2.000\n",
       "Prob(Omnibus):                  0.610   Jarque-Bera (JB):                0.999\n",
       "Skew:                           0.006   Prob(JB):                        0.607\n",
       "Kurtosis:                       2.991   Cond. No.                         1.05\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = smf.ols('X ~ Y', model).fit()\n",
    "ols.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.60281\n",
       "dtype: float64"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols.predict({'Y': 2})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use a similar approach to estimate the values of variables in the middle of a causal chain when we know the values before and after. For example, if we want to estimate the value of Y given that we observe specific values for X and Z, all we have to do is fit a linear model for Y as a function of X and Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>Y</td>        <th>  R-squared:         </th>  <td>   0.105</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   0.105</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>   5836.</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sat, 25 Jul 2020</td> <th>  Prob (F-statistic):</th>   <td>  0.00</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>17:37:39</td>     <th>  Log-Likelihood:    </th> <td>-1.4135e+05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>100000</td>      <th>  AIC:               </th>  <td>2.827e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td> 99997</td>      <th>  BIC:               </th>  <td>2.827e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -0.0026</td> <td>    0.003</td> <td>   -0.815</td> <td> 0.415</td> <td>   -0.009</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X</th>         <td>    0.3317</td> <td>    0.003</td> <td>  105.493</td> <td> 0.000</td> <td>    0.326</td> <td>    0.338</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Z</th>         <td>    0.0656</td> <td>    0.003</td> <td>   20.906</td> <td> 0.000</td> <td>    0.059</td> <td>    0.072</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 2.830</td> <th>  Durbin-Watson:     </th> <td>   2.002</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.243</td> <th>  Jarque-Bera (JB):  </th> <td>   2.821</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.013</td> <th>  Prob(JB):          </th> <td>   0.244</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.007</td> <th>  Cond. No.          </th> <td>    1.02</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      Y   R-squared:                       0.105\n",
       "Model:                            OLS   Adj. R-squared:                  0.105\n",
       "Method:                 Least Squares   F-statistic:                     5836.\n",
       "Date:                Sat, 25 Jul 2020   Prob (F-statistic):               0.00\n",
       "Time:                        17:37:39   Log-Likelihood:            -1.4135e+05\n",
       "No. Observations:              100000   AIC:                         2.827e+05\n",
       "Df Residuals:                   99997   BIC:                         2.827e+05\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.0026      0.003     -0.815      0.415      -0.009       0.004\n",
       "X              0.3317      0.003    105.493      0.000       0.326       0.338\n",
       "Z              0.0656      0.003     20.906      0.000       0.059       0.072\n",
       "==============================================================================\n",
       "Omnibus:                        2.830   Durbin-Watson:                   2.002\n",
       "Prob(Omnibus):                  0.243   Jarque-Bera (JB):                2.821\n",
       "Skew:                           0.013   Prob(JB):                        0.244\n",
       "Kurtosis:                       3.007   Cond. No.                         1.02\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = smf.ols('Y ~ X + Z', model).fit()\n",
    "ols.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now that we have a description of Y as a function X and Z, we can simply plugin the numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.525892\n",
       "dtype: float64"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols.predict({'X':1, 'Z':3})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Blocking variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An interesting point to make clear (and that we will delve into more deeply in Chapt 2) is that to determine the value of a given variable, we only need the know the values of it's parents. \n",
    "\n",
    "Let's say we observe that Y=3 and are interested in finding out the value of Z. We can simply do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.197132\n",
       "dtype: float64"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = smf.ols('Z ~ Y', model).fit()\n",
    "ols.predict({'Y': 3})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, let's say that we also observe that X=1. In this case, we would do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>Z</td>        <th>  R-squared:         </th>  <td>   0.005</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   0.005</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>   244.6</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sat, 25 Jul 2020</td> <th>  Prob (F-statistic):</th>  <td>1.07e-106</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>17:37:39</td>     <th>  Log-Likelihood:    </th> <td>-1.4195e+05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>100000</td>      <th>  AIC:               </th>  <td>2.839e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td> 99997</td>      <th>  BIC:               </th>  <td>2.839e+05</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -0.0026</td> <td>    0.003</td> <td>   -0.823</td> <td> 0.410</td> <td>   -0.009</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Y</th>         <td>    0.0664</td> <td>    0.003</td> <td>   20.906</td> <td> 0.000</td> <td>    0.060</td> <td>    0.073</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X</th>         <td>    0.0007</td> <td>    0.003</td> <td>    0.217</td> <td> 0.828</td> <td>   -0.006</td> <td>    0.007</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.440</td> <th>  Durbin-Watson:     </th> <td>   2.008</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.487</td> <th>  Jarque-Bera (JB):  </th> <td>   1.444</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.003</td> <th>  Prob(JB):          </th> <td>   0.486</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.983</td> <th>  Cond. No.          </th> <td>    1.39</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      Z   R-squared:                       0.005\n",
       "Model:                            OLS   Adj. R-squared:                  0.005\n",
       "Method:                 Least Squares   F-statistic:                     244.6\n",
       "Date:                Sat, 25 Jul 2020   Prob (F-statistic):          1.07e-106\n",
       "Time:                        17:37:39   Log-Likelihood:            -1.4195e+05\n",
       "No. Observations:              100000   AIC:                         2.839e+05\n",
       "Df Residuals:                   99997   BIC:                         2.839e+05\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.0026      0.003     -0.823      0.410      -0.009       0.004\n",
       "Y              0.0664      0.003     20.906      0.000       0.060       0.073\n",
       "X              0.0007      0.003      0.217      0.828      -0.006       0.007\n",
       "==============================================================================\n",
       "Omnibus:                        1.440   Durbin-Watson:                   2.008\n",
       "Prob(Omnibus):                  0.487   Jarque-Bera (JB):                1.444\n",
       "Skew:                          -0.003   Prob(JB):                        0.486\n",
       "Kurtosis:                       2.983   Cond. No.                         1.39\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = smf.ols('Z ~ Y + X', model).fit()\n",
    "ols.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.197201\n",
       "dtype: float64"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols.predict({'Y': 3, 'X':1})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which returns the same result (up to numerical precision). While counterintuitive, this result is simple to understand. In our causal model, the value of Z is determined only by the value of Y. Once Y is determined, the factors that might have led to its value (a statistical fluctuation in Uy, or the value of X) are irrelevant. All that Z \"sees\" is the current value of Y and Y *shields* Z from everything that happened before in its causal chain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this simple observation lies much of the power of this approach!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Product Decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = CausalModel('dags/Primer.Fig.1.10.dot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAYAAACmKP9/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAuIwAALiMBeKU/dgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxddZ3/8fe5N1vTLSlp7UKXJGVpoNCN1R+0iEgf6E+BAqMoMlIEygiio6ODitQNHETZZBVwHEd2ERzGuiCgbKE0tLK1tGkLaRsoIU3TLdu9Z/4It/Q0SbOdc77fc87r+Xj48MG3zb3f3N6cT17nbo7ruq4AAECspExvAAAA+I8BDwBADDHgAQCIIQY8AAAxxIAHACCGGPAAAMQQAx4AgBhiwAMAEEMMeAAAYogBDwBADDHgAQCIIQY8AAAxxIAHACCGGPAAAMQQAx4AgBhiwAMAEEMMeAAAYogBDwBADDHgAQCIIQY8AAAxxIAHACCGGPAAAMQQAx4AgBhiwAMAEEMMeAAAYogBDwBADDHgAQCIIQY8AAAxxIAHACCGGPAAAMRQnukNAPCX67pq2CnVbnG1Yaurtoyr1ozUnun88/y0VJiWCtKO9h/pqLLUUVmx5DiO2Y0D8BUDHoi4hh2u1mxxtbYxq9otrmobXTW39u8yRhRKlaM6h33FqJSmljoqG8rAB6LMcV3XNb0JAP3TnnH1bF1WS9ZktaohmB/hg8oczZ+a0rETU8pPM+yBqGHAAxGyeburP9Vm9Je12X5X+kCNKJQ+WpHSxyrTGjOMQQ9EBQMeiICX38nq96syWrbJlakfWEfS7PGO/v9BaU3/EM/PBWzHgAcstqPN1V0vZfTEuqzprXicUJ7SeTPTGlpA0QO2YsADllq2KatblnaocZfpnXRv1BBp0RF5mj2emgdsxIAHLGNrtfeEmgfsxIAHLLK8PqubXrC32nsyaoj0pSPzNGMcNQ/YggEPWOKvazP6+QsZY0+iGyxH0peOSuuE8rTprQAQAx6wwv+syuiulzKmt+GLhbPS+viBDHnANAY8YNjDr2f0XyviMdxzzjk8rdOmMeQBk3jADDDoj2viN9wl6b9WZPTHNfH7voAoYcADhvz9zYxufzG+Q/D2FzN6+s34fn+A7RjwgAF1W13dWB3dJ9T1hSvphuqM6rbG+bsE7MWAB0KWybq6sbpDHdF4mfugdGSlm6o7lMky5IGw8XGxQMgeWZXVmsaBD7wd767Xny49yMcd9e7gBd/WtDO+M6CvXd3o6tFVWZ50B4SMggdCVLfV1b0vJ+9x6XtfzmgDp+qBUDHggZBksq5uSsip+b21Z6UbOVUPhIpT9EBIHlmV1epBnJrPSaULNHLyYT7sqKsd76xVR8v2LuvFZZMGfdmcqgfCxRvdACHY0ebqi4+2q6XD9E569t6qZ/X3739Mbqbdsz7hmDN15KW/9uU6ivKkOz6ZzwfTACHgFD0QgifXZ60e7rsaN6r6Z5/uMtxHTJquWRfe7tv1tHRIT61P4GMUgAEMeCBgrutqicXv6pZpb1X1T/9JrVvf8aznDxulo//1AeUVFvt6fUvWZMSJQyB4DHggYK9sdrWx2fQuerb8zku0pXapZ81JpXXkpb/W0DHlvl/fhmbp1c0MeCBoDHggYDbXe+2Sn+utp/6zy/ohn/mBxkw/MbDrXbKG0/RA0BjwQIAad7mq3mBnrTa89je9/Ot/67K+/7H/pAM+8dVAr7t6Q1aNu+y8XYC4YMADAfpLbVY2vvR753t1qr7+bLkZ7zP/Rk4+XDMvuDXw68+40uNrqXggSAx4IECvbLZviGXaWlR97Vlqa37Xs14wbD8d9a/3+/6kup688o59tw0QJwx4ICCu62rtFvvy/aU7FqlpXY1nzUmldcSXf62ho6eEto+1W1yeTQ8EiAEPBOTt7dLO9t7/XphWP3ad6p7+TZf1Q87+kcYc+pFQ97KjXXqn65vmAfAJAx4ISG2jXaegN7/yV736m8u7rE/88Gd0wMcvM7Ajac0Wu24jIE4Y8EBAai06Pb/j3fVaev3n5Ga9L9kbOWWGZl5wi6FdSWt9eG9+AN1jwAMBqbVkeHW07lT1tWeqbft7nvWC4WU6+qv3K10wxNDO7LmNgDhiwAMBsOkJdjW3XaCtb/7Ds+ak83Tkl/9bxaMnG9pVJ55oBwSHAQ8EoCNrxxPs3nj0J9r43ANd1g89+yqNPmRe+Bvay472ztsKgP8Y8EAA2ix4d9p3/vFnvXbfFV3WJx73WU095VIDO+qeDbcVEEcMeCAApofW9ndqtfSGc7o8qa6kfJZmnn+zoV11r52CBwLBgAcCYHJodbRsV/W1Z6p9xxbPesGI0Trqq/cpXVBkaGfda6fggUAw4IEA5Bn8yVp2y/lqrnvVs7b7SXVlkwztqmcmbysgzvjRAgJQmDZzvasevlqbXni4y/r0z/1Yo6vmGthR7woM3VZA3DHggQDkGxhab7/0B732wOIu65OOP0eV878U/ob6iAEPBIMBDwQgPyUNyQvv+rbVv6EXbzpXcr0P/pdUzNaMhTeFt5F+GpLHKXogKPxoAQFwHEflpU4o19W+a1vnk+p2bvWsF44YY+WT6vZUUerIccK5nYCkYcADAakcFfzgcl1Xy27+grZtXOlZd9J5OvKy36h4v4mB72EwKkK4jYCkYsADAakIoeBXPvQD1b/4+y7r08+5RmXTjgv8+gerMqSzHEASMeCBgEwdFeyPV/2Lv9fK3/6wy/qkuZ9X5ckXB3rdfgn6NgKSjJ8uICDjhgf3RLttG1fqxZvPk/b6oJbSyiOsflLdnobkSWOHm94FEF8MeCAgqYCeaNe+c6uev/YMdexq9qwXjvyQjvrKvUrnF/p+nUGoKHWU4gl2QGBCfCEPkDyHjHH02rv+fhzqil9+VdvrV3dZd92snrvmNF+va08fuXqpr5d3yBiGOxAkBjwQoBMr0nrw1az8HPG7Gt7qdr2t+V21Nb/r4zUFx1HnbQMgOJyiBwI0ZqijORMo1b0dMcHR6KHcLkCQGPBAwOZPpVT3xm0CBM9xXdffBwgBeGRdV196rF1vbze9EzuMGybd+PF8nmAHBIyCBwKWchydTLHudvLUNMMdCAEDHgjBR8pTRj5hzjb5aemEcg47QBj4SQNCMLzQ0WkH8+N22sEpDS+k3oEwcMQBQrKgKq0pJckdblNKHC2o4jQGEBYGPBCS/LSjS45KK53AGZ92pEuOSis/id88YAgDHghReWlKC6qS92O3oCql8tLkfd+ASfzEASFL2ql6Ts0DZjDggZAl6VQ9p+YBcxjwgAHlpSldMCf+VXvBnDSn5gFD+MkDDDmpMq3PHx7fIf/5GWmdVBnf7w+wHQMeMOjUaWmdEcMn3Z1RldKpBzPcAZN4L3rAMNd19cCrWd37Ssb0Vnzxmemdv7Q4vB0tYBQDHrDEY29kdGdNtIf8+bPSOuVAyh2wAQMesMhzdVndurRD29pM76R/hhdIi47I09ET4/dwAxBVDHjAMk0trm57sUPVG6Lxo3n0/o4umJOnkiJOyQM2YcADFnJdV8+8ldUdyzLa1tb537Y8pp3by/AC6Ytz0vrwRB5vB2zEgAcs1tTi6ruP1qkuO1aS2UG/53VPSr2tKz85kWoHLMaAByyWzWY1a9Ysbc6frKmnXKKyaccb3U/D63/Tmv+9UR/qeEs1NTWUO2AxBjxgsYcfflinn3767v8evn+VLrzmQW3Mq1RLRzh7KMqTJrTX6rZ/O0PbNrzm2dupp54aziYA9BsDHrBUrt5XrFixe23atGl6+eWX1ZZN6an1Wf1hTVZ1W4P5EZ400tH8qSnNnZJSQSqrQw89VCtXrtz95zNmzKDiAYsx4AFL7V3vknTPPffo05/+9O7/dl1Xdc2uahs/+N+6Jldt/Xw5fUFaKi9xVDnqg/9NHOF4hvc999yjs88+u8seqXjATgx4wEL7qvd0et9vJJPJutrQ7Grtls7/b+2Q2jLaPfQL0p3/K8yT9h/hqKLU0f4jHKVT+y7xTCZDxQMRwoAHLNSXejeBigeigwEPWGYw9R40Kh6IDt5XErDMI4884hnuknTFFVcYH+6SlE6ndcUVV3jWli9frkceecTQjgD0hIIHLGJzvedQ8UA0UPCARWyu9xwqHogGCh6wRBTqPYeKB+xHwQOWiEK951DxgP0oeMACUar3HCoesBsFD1ggSvWeQ8UDdqPgAcOiWO85VDxgLwoeMCyK9Z5DxQP2ouABg6Jc7zlUPGAnCh4wKMr1nkPFA3ai4AFD4lDvOVQ8YB8KHjAkDvWeQ8UD9qHgAQPiVO85VDxgFwoeMCBO9Z5DxQN2oeCBkMWx3nOoeMAeFDwQsjjWew4VD9iDggdCFOd6z6HiATtQ8ECI4lzvOVQ8YAcKHghJEuo9h4oHzKPggZAkod5zqHjAPAoeCEGS6j2HigfMouCBECSp3nOoeMAsCh4IWBLrPYeKB8yh4IGAJbHec6h4wBwKHghQkus9h4oHzKDggQAlud5zqHjADAoeCAj1/gEqHggfBQ8EhHr/ABUPhI+CBwJAvXdFxQPhouCBAFDvXVHxQLgoeMBn1HvPqHggPBQ84DPqvWdUPBAeCh7wEfXeOyoeCAcFD/iIeu8dFQ+Eg4IHfEK99x0VDwSPggd8Qr33HRUPBI+CB3xAvfcfFQ8Ei4IHfEC99x8VDwSLggcGiXofOCoeCA4FDwwS9T5wVDwQHAoeGATqffCoeCAYFDwwCNT74FHxQDAoeGCAqHf/UPGA/yh4YICod/9Q8YD/KHhgAKh3/1HxgL8oeGAAqHf/UfGAvyh4oJ+o9+BQ8YB/KHign6j34FDxgH8oeKAfqPfgUfGAPyh4oB+o9+BR8YA/KHigj6j38FDxwOBR8EAfUe/hoeKBwaPggT6g3sNHxQODQ8EDfUC9h4+KBwaHggd6Qb2bQ8UDA0fBA72g3s2h4oGBo+CBfaDezaPigYGh4IF9oN7No+KBgaHggR5Q7/ag4oH+o+CBHlDv9qDigf6j4IFuUO/2oeKB/qHggW5Q7/ah4oH+oeCBvVDv9qLigb6j4IG9UO/2ouKBvqPggT1Q7/aj4oG+oeCBPVDv9qPigb6h4IH3Ue/RQcUDvaPggfdR79FBxQO9o+ABUe9RRMUD+0bBA6Leo4iKB/aNgkfiUe/RRcUDPaPgkXjUe3RR8UDPKHgkGvUefVQ80D0KHolGvUcfFQ90j4JHYlHv8UHFA11R8Egs6j0+qHigKwoeiUS9xw8VD3hR8Egk6j1+qHjAi4JH4lDv8UXFAx+g4JE41Ht8UfHAByh4JAr1Hn9UPNCJgkeiUO/xR8UDnSh4JAb1nhxUPEDBI0Go9+Sg4gEKHglBvScPFY+ko+CRCNR78lDxSDoKHrFHvScXFY8ko+ARe9R7clHxSDIKHrFGvYOKR1JR8Ig16h1UPJKKgkdsUe/IoeKRRBQ8Yot6Rw4VjySi4BFL1Dv2RsUjaSh4xBL1jr1R8UgaCh6xQ72jJ1Q8koSCR+xQ7+gJFY8koeARK9Q7ekPFIykoeMQK9Y7eUPFICgoesUG9o6+oeCQBBY/YoN7RV1Q8koCCRyxQ7+gvKh5xR8EjFqh39BcVj7ij4BF51DsGiopHnFHwiDzqHQNFxSPOKHhEGvWOwaLiEVcUPCKNesdgUfGIKwoekUW9wy9UPOKIgkdkUe/wCxWPOKLgEUnUO/xGxSNuKHhEEvUOv1HxiBsKHpFDvSMoVDzihIJH5FDvCAoVjzih4BEp1DuCRsUjLih4RAr1jqBR8YgLCh6RQb0jLFQ84oCCR2RQ7wgLFY84oOARCdQ7wkbFI+ooeEQC9Y6wUfGIOgoe1qPeYQoVjyij4GE96h2mUPGIMgoeVqPeYRoVj6ii4GE16h2mUfGIKgoe1qLeYQsqHlFEwcNa1DtsQcUjiih4WIl6h226q/jDDz9cNTU1SqVoJdiHeyWsRL3DNt1V/IoVK/Too48a2hGwbxQ8rEO9w1Y8Fo8ooeBhHeodtuKxeEQJBQ+rUO+wHRWPqKDgYRXqHbaj4hEVFDysQb0jKqh4RAEFD2tQ74gKKh5RQMHDCtQ7ooaKh+0oeFiBekfUUPGwHQUP46h3RBUVD5tR8DCOekdUUfGwGQUPo6h3RB0VD1tR8DCKekfUUfGwFQUPY6h3xAUVDxtR8DCGekdcUPGwEQUPI6h3xA0VD9tQ8DCCekfcUPGwDQWP0FHviCsqHjah4BE66h1xRcXDJhQ8QkW9I+6oeNiCgkeoqHfEHRUPW1DwCA31jqSg4mEDCh6hod6RFFQ8bEDBIxTUO5KGiodpFDxCQb0jaah4mEbBI3DUO5KKiodJFDwCR70jqah4mETBI1DUO5KOiocpFDwCRb0j6ah4mELBIzDUO9CJiocJFDwCQ70Dnah4mEDBIxDUO+BFxSNsFDwCQb0DXlQ8wkbBw3fUO9A9Kh5houDhO+od6B4VjzBR8PAV9Q7sGxWPsFDw8BX1DuwbFY+wUPDwDfUO9A0VjzBQ8PAN9Q70DRWPMFDw8AX1DvQPFY+gUfDwBfUO9A8Vj6BR8Bg06h0YGCoeQaLgMWjUOzAwVDyCRMFjUKh3YHCoeASFgsegUO/A4FDxCAoFjwGj3gF/UPEIAgWPAaPeAX9Q8QgCBY8Bod4Bf1Hx8BsFjwGh3gF/UfHwGwWPfqPegWBQ8fATBY9+o96BYFDx8BMFj36h3oFgUfHwCwWPfqHegWBR8fALBY8+o96BcFDx8AMFjz6j3oFwUPHwAwWPPqHegXBR8RgsCh59Qr0D4aLiMVgUPHpFvQNmUPEYDAoevaLeATOoeAwGBY99ot4Bs6h4DBQFj32i3gGzqHgMFAWPHlHvgB2oeAwEBY8eUe+AHah4DAQFj25R74BdqHj0FwWPblHvgF2oePQXBY8uqHfATlQ8+oOCRxfUO2AnKh79QcHDg3oH7EbFo68oeHhQ74DdqHj0FQWP3ah3IBqoePQFBY/dqHcgGqh49AUFD0nUOxA1VDx6Q8FDEvUORA0Vj95Q8KDegYii4rEvFDyodyCiqHjsCwWfcNQ7EG1UPHpCwScc9Q5EGxWPnlDwCUa9A/FAxaM7FHyCUe9APFDx6A4Fn1DUOxAvVDz2RsEnFPUOxAsVj71R8AlEvQPxRMVjTxR8AlHvQDxR8dgTBZ8w1DsQb1Q8cij4hKHegXij4pFDwScI9Q4kAxUPiYJPFOodSAYqHhIFnxjUO5AsVDwo+ISg3oFkoeJBwScA9Q4kExWfbBR8AlDvQDJR8clGwccc9Q4kGxWfXBR8zFHvQLJR8clFwccY9Q5AouKTioKPMeodgETFJxUFH1PUO4A9UfHJQ8HHFPUOYE9UfPJQ8DFEvQPoDhWfLBR8DFHvALpDxScLBR8z1DuAfaHik4OCjxnqHcC+UPHJQcHHCPUOoC+o+GSg4GOEegfQF1R8MlDwMUG9A+gPKj7+KPiYoN4B9AcVH38UfAxQ7wAGgoqPNwo+Bqh3AANBxccbBR9x1DuAwaDi44uCjzjqHcBgUPHxRcFHGPUOwA9UfDxR8BFGvQPwAxUfTxR8RFHvAPxExccPBR9R1DsAP1Hx8UPBRxD1DiAIVHy8UPARRL0DCAIVHy8UfMRQ7wCCRMXHBwUfMdQ7gCBR8fFBwUcI9Q4gDFR8PFDwEUK9AwgDFR8PFHxEUO8AwkTFRx8FHxHUO4AwUfHRR8FHAPUOwAQqPtoo+Aig3gGYQMVHGwVvOeodgElUfHRR8Jaj3gGYRMVHFwVvMeodgA2o+Gii4C1GvQOwARUfTRS8pah3ADah4qOHgrcU9Q7AJlR89FDwFqLeAdiIio8WCt5C1DsAG1Hx0ULBW4Z6B2AzKj46KHjLUO8AbEbFRwcFbxHqHUAUUPHRQMFbhHoHEAVUfDRQ8Jag3gFECRVvPwreEtQ7gCih4u1HwVuAegcQRVS83Sh4C1DvAKKIircbBW8Y9Q4gyqh4e1HwhlHvAKKMircXBW8Q9Q4gDqh4O1HwBlHvAOKAircTBW8I9Q4gTqh4+1DwAXj++ee1cOFCLV68WPX19d3+HeodQJz0t+Lr6+u1ePFiLVy4UM8//3wYW0wcCt5n9fX1mjp1qnbu3ClJKioq0oUXXqhvfOMbGjdunCTqHUA89aXi6+vr9eMf/1i33XabWlpaJEnFxcVas2bN7mMk/EHB+6y6unr3cJeklpYWXX/99aqoqNBll12m+vp66h1ALO2r4uvr63XZZZepoqJC119//e7hLkk7d+5UdXV12NuNPQreZw888IDOOuusHv+8qKhIw4YNU0NDw+416h1AXHRX8WVlZdq+fbtnqO/t/vvv15lnnhnGFhODgg9ZS0uLZ7hL1DuA+Oiu4hsaGvY53BEMBrxhjuPoueee6/HJeAAQJfX19Xr++ed55rwFGPCGua6rG264wfMYPQBEzZ6Psd9www3i0V/zGPCW6O7JeABgu309eQ5mMeAtkxv0VVVVWrp0qentAECPli5dqmnTpjHYLcWAt1RTU5Ouu+4609sAgB5dd9112rp1q+ltoAcMeIuVlZWZ3gIA9IhjlN0Y8JaaM2eOFi9ebHobANCjK6+8UrNnzza9DfSAAW+hSy+9VE8//bRKSkpMbwUAelRaWqpnnnlGl1xyiemtoBsMeIuMHDlSDz30kK6//noVFhaa3g4A9KqwsFA33HCDHnzwQY0YMcL0drAHBrwl5syZo5qaGp1++ummtwIA/bZgwQK99NJLnLK3CAPeArlT8hUVFaa3AgADVlFRwSl7izDgDeKUPIC44ZS9PRjwhnBKHkCcccrePAa8AZySB5AEnLI3iwHvs2w22+OfcUoeQNL09ZQ9H07jvzzTGzBlW6ur2i2u1m5x1bTLVXtWauuQOrJSXkrKT0sFaWlEoaOKUY4qSx2NLOr94w83bNjQ7fqcOXN03333Ue0AEmnBggWaOXOmzjrrLC1btqzLn9fV1fXpcra2vH/sbnTV3OqqLSO1Zz44dhfkSfkpqWSIo4rSzmP38MJkfnRtIgb8znZXb7zXeYeo3ZJVbaOrzTv6fzllxVLlKEeVpSlVjHJ04H6OhhV47zhjx47t8nWXXHKJrrnmGqodQKLlTtl/7Wtf00033eT5s+6Ondvbuh67G3b2/3rHDO167C7Oj//Qd9wYnxepbcxqyZqs/v5mVm0Z/y8/PyUdOyml+VNTOnA/R47jKJvNasaMGXr55ZeVTqf1k5/8RJdddpn/Vw4AEfazn/1MX//615XJZDR9+nQtX75cqVRKrts51JesyeqZt7Lq6PlRzwErSEvHTe48dleOiu8j1bEb8G0ZV8+81TnYV78X3rdWXupo/tSUjpucUlGeo8bGRpWUlCiViu+dBwAGI5vNqqmpSaNGjVJLh6u/v9l57F63Jbxj9wH7dR67PzwppYJ0vKo+NgN+a4urR1Zm9PjarLa1mdtHcb50QnlKp01La9SQeN1ZAMBvjbtcPfx6Rk+sy2pnu7l9DC+QTqxI6VMHp/v0fKsoiMWAf7Yuq9tf7FBzq+mdfGBovrRwVlpzp6TkOPG4swCAX1zX1ZPrs7qrJqMdBgf73kYUShfOydMxE6N/9jXSA35ri6s7lmX0bF0AD9L4ZM54RxcdkUfNA8D7Gne5unVph17cZO/4OXZiSl+cHe2aj+yAt7Hae0LNA4C91d6TqNd85Ab8jjZXtyy1u9p7Mme8oy8dlacRCX1NJoDkam51dVO13dXek2MnprToiLSGFkTr2B2pAd/U4ur7T3ZoXVNkttzFhBHSd+flq6w4WncUABiohp2uFj/Zro3NpncycOWljr4zN08lETplH5kB/+6OzjvIpm2mdzJ4o4ulK0/I17jh0bmjAMBAbNrmavET7Xp3AG9QY5sJw98PtKHROHZHYsA37HD1rcfjcQfJKSmSfnBivsYz5AHE1KZtrr79eLuaWkzvxD+ji6UfnhiNIW/9MweaWjrLPU7DXZKaWqTFT7SrYaf1v18BQL817HR15RPxGu6S9O5OafGT7draYv+x2+oBv6PN1fef6tDGGJyW707ujtLcav8dBQD6qrm1M8wG8r7xUbBxm/S9pzq0o83uY7fVA/6OZZlQ37LQhI3N0k3VHXxUIoBYcN3OZ8tH+Ql1fbFui6tfLAvgQ058ZO2nyb2wIau/vTm4l8K9/sD3tPK3P/SslR5wlOZe+YScVHrAl+u6rp696uPa/PLjnvXK+f+iw879ab8v78VNrp5an9W88oHvCQBs8OT67IBfCudmM3rqirnaUrvUs37wgu9o2hnf9mN7am1+V49/fZZamzd71o/++m81btbH+3VZT72Z1TGTsjpygp2tbOWutrW6uvXFjkFfzkGn/btGTj7cs7ZldbVW/8/PBnW56/58W5fhPmzcATrkMz/s4St6d2dNRo27qHgA0dW4y9VdNQOvWieV1uyL71S6YIhnfdXvrlbTupcGuz1J0vI7L+ky3CfP++d+D/ecW5d2aJulD7NaOeDvrMn48sSMVF6+Zi+6U6m8As/66w9+T811rw3oMre/U6tXfnO5Z81JpTV7Udc7ZX/saO+8o3CqHkAUuW7n288O9h3qho8/SFWf/r73sjPtevHm85RpH9xbl9Y9fY82vfCwZ624bLKmn3PNgC+zqUWD+qUmSNYNeD9Oze9p5OTpOniB99ROtr1Vy25ZqGymf2cJ3GxWy245X5nWHZ71qZ/4ikYdcNSg95o7VQ8AUTOYU/N7q5z/JZVVzfWsbdvwml67/7sDvsxdjZu04pdf8S46jmZddLvyi0cM+HKlzlP1L2y079ht1YD369T83g785NdUWnmEZ61pXY3e+N2P+3U5qx/7mRpXPetZGzHxEFWdOfA73d44VQ8gagZ7an5vjuNo9kV3KG/IcM/6mseuV8PKZwZ0mS/dcZHad2zxrFWe/C8afci8gW7Tw8ZT9VYN+Mfe8OfU/N6cVFqzFt2hVH6RZ33lw1epaf2KPl1Gc91rev2Bxd7LTXf/EMBg7GiXHn7dztM9ANCdh1/3/8NjikdP1vRzfuJddLOquWWhOlq29+uy1j1+p95Z/kfP2rDxB+qQz/xgsNvcramlc4bZxJoB35F19S23qaYAAA+ASURBVKfa4E5xjJgwTVVnXelZczPtWnbzecp2tO3za7OZDi275Txl93r856BTv6mS8pl+b1VPrMuqpcOu3wQBoDstHa6eWBfMsXvKCf+ssbNO8azt2LxOL//6G32+jB3vrtcr/+39+53Pm7prUM+b6s6fa7PqyNpz7LZmwL+wwQ38HY+mnvJljTroWM9ac90rev3B7/fwFZ1WPXxVl2dwlpTP1EGnfdP3PUrSznbpaR+fhwAAQfn7m1ntDPCjX2d+8RYVDNvPs7b+8V90KfLuuK6rmlsvUMcu77ulHfipf9OoqUf08FUDt6VFemEjA76LJWuCP7XhpFKafdEdShcWe9ZX//5aNa5Z2u3XbFlbo1W/u9qzlsov1OxFdymVDu5tBP6wJssz6gFYzXVdLVkdbIwUlYzVjIU3dlmvuf0itW3f0s1XfKB2yU1qeO0pz9rIyYfr4NO/5ese97RktT2n6a0Y8HVbXb2yOZxhNmzs1C6vV3ezGS27ZaEybd5TCJn2Vi27+Ty5ez3bftoZV2jExKpA97lui6vVjQx4APZ64z03lI/vnnD0Au1/7D951lq2bNKKX17W49dsq39Dr937Hc9aKr9Qc/7lbqXy8gPZpyS9stnVhq12HLutGPBh1PueKj62SGVV8zxr2zet0mv3XeFZe+3+72rbxtc9a6MOOFoHfGKvl1oEJOjfjAFgMJasCe8YdfgXrldR6XjP2oZn7tXG6t92+bud0Xa+Mm27POudcXZIoPuUwp9pPTE+4Fs7XD0Z0BM0etL5Eozbu74E4w837n4JRsPKZ7Tmses9f54uLNbsRb8Y1Nvc9sczb2W13fIPMwCQTNvbXD3zVnjH7oJhpZp5wS1d1pffeYlamt7xrL3x+59qy+pqz9qog47VAZ/4aqB7zHlyfVatFjxR2viAX9/kapf/L33vVfHoyTr0s3u9Dt7NqubW89W6dbNqbj1fcr133kM+/QMNG3dAaHtsz0qr3zN/JwEQTdXV1Zo9e7aqqqr0ox/9SNu2+ffRnG+856oj5JOMY2fM15QTF3rW2rY16KU7Fu3+7+a6V7Xywe95/k66cOj7cRbOyNvZ3jnbTDM+4GsNPs5cfuJCjTnsJM/ajnfW6vFvHqEd76z1rJdVzVPFyReHuT1JZm8fANF2zjnnqKamRq+//rq+9a1vacqUKbrqqqt8GfRrDR2bpn/uP1Q8Zopn7e2ax/Tmk/+pbEfnW9ru/dLnQz97lYZ9qDLEXUprLfgkVPMD3vCNMOvC25Q/tMSz1tr0tue/84YM1+yLbpfjOGFuTZK0dguPwwPov/b2dq1evdqz1tjYqMsvv9yXQV9r6NiUVzRMsy/6heR4x9c/fvU1rbj7Um1dv9yzPuawj6ripAvD3KIkO+LM/IA3fCMMGTVBh33+2n3+nUM/+2MVj54c0o681lhwJwEQL34MepPH7rJpx2nqKZd61jp2NWv9X+/yrOUPLdGsC24Lc2u7mZ5tkuEB39rhakOz+Rth0vGf05jDPtrtn31oxskq3+sxnzA17JSaLXt/YwDxMNBBv7XFVcPOgDfXi6qzFmv4hGn7/DuHnfszDdlv/5B25FXX7Bp/op3RAb++yZUN7+rX2tygrW++3O2fjT/y9JB305UNvwkCiK/+DnrTD61KUrqgSLMvvktOD284Nv7IUzXpuLND3tUHsq75J9oZHfA2PAlBklbcdalat77T7Z+9eu+31dr8bsg78rLldgIQb30d9KaeYLe30opZOujUrm8ZXjC8TDMW3mRgR16mj91GB/zWFvN3kg3P3q+N1Q/1+Odtze9qxV1fDnFHXXGKHkCYehv0Nh2TyqqO77I2Yv8qFY4YbWA3Xs2tvf+dIAX3Zup90Gb4zX5amt7W8ru7Du/84pFq37l1939vrH5IG59/SBOOXhDm9nb741+e1NM33WPkugFEUzY7+Ge55wb94sWLdfrpp+tXv/qV8vLyjB+7o6ItY/YXIaMDvt3wK8BeuuNitW9v9KyVf/SLGnPYSar+6Vme9eV3f1llVccb+a1wde061dx1V+9/EQAC0NraqnvuuUeNjY1asmSJ2hnwfWL6FyGjp+jDfhekPb351K/0ds1jnrXiMVN06Gev1vgjPtXlgw3amt/V8ru8L8sISyq/0Mj1AsCeli7t/NRNk8fuKDF9Oxkd8HmGrn3Xexv08q++5l10HM2+8A7lFQ2TJB3+hetUWDLW81c2Vf9WG557MKxt7pZtN/xADgBIOv74zse7TR27o8b07WT0FH1BOJ/Z0kXN7Rd6HmOXpMqTL/Y8WaNg2CjNWHijqq890/P3Vtz9ZY2uOl6FI8eEsldJOnTagZr7lXA+wQ5APGQyGd1www2+XNawYcN07rnn7r68AqOTIzoKDc24HKP/TPkGfrtZ9+fbtfkff/GsDR07VVWf/kGXvzt+zic18f+drbqnf7N7rW1bg5bfdamO+sq9ge81Z97xH9YXZs4N7foARF97e/ugB/yECRN0+eWXa+HChSos/OChQhPH7ijKT4f/9uZ7MvrPVDok3G9+x+Z1euU3/+5ddFKavegXyiss7vZrDjv3pyoqGedZ2/TCw9rw3ANBbbOLkUVm7yQAkmXChAn6+c9/rtraWl188cWe4S5JJSEfu6OqpMjs9Rsd8BWl4d1JXNdVza0XqKNlu2d96se/rP0OPKbHrysYVqoZ53d9w4QVd1+m1q2bfd9ndypDvJ0AJFdvgz0nzGN3lJm+nYwO+MkljsI6g1H7hxvV8PrfPGvDJxysqjOv7PVrx83+hCYe9znPWuep+kv83GKPTN9JAMRbXwd7DtHRu7TTOeNMMjrgC9KOJo4M/gbYVv+GXrvvCs+ak0pr9qI7lS7o2zmUw869VkWl4z1rm174nTY8e79v++zOh4ZKwwv5YQLgv/4O9pzhhY7GDA14cxE3aaSjgiQ/Bi9JlaOCvQHcbEbLbjlfmbZdnvUDPvk1lVbO6fPlFAwt0czzf95lfcXdl6mlqfv3sfdD5Sjj/0QAYmagg31PQR+7o86G28f49Aj6VM8bv/+ptqyu9qyNmHiopi34dr8va+ysUzRp7uc9a23b3wv0VH2FBXcSANGTn5+vqqoqz5ofgz2nstT4+LAaA17B3gjNda9q5YPf86w56XzNXnSnUnkFA7rM6edco6JREzxr9UsfUd0z9w14n/vCY10ABuqee+7RCSecoGOOOca3wZ5DfOybDc+dMv52BVNKHA0rkLa3+Xu52Y52LbtlobId3gs+6NRvqqR8xoAvt2BoiWZ+8WY99+NPedb/8cuvaPQh81RU8qEBX3aX60pLB+xn/k4CIJoOO+ww/fWvfw3ksg/cz1FB2vz7rdtoWEHnbDPNeMHnpx2dUO7/Nlb97mo1rXvJs1ZSPlMHndb1s4P7a+yM+Zo87589a23b39PyO7806Mve0/GTUyrON38nAYC9Fec7Om6y8RFipY+Up4y/yY1kwYCXpJOn+vt+fk3rXtKq313tWUvlFXSemk/7c9Ji+jn/oSGj9ves1b/4qOqe8e8d7k6easU/DwB0az7HqG75PdMGynFd1+wH1r7ve0+2a/nbVmzFCgfu5+jqk/JNbwMA9ukbf27X6vc4dufMGOvoinl2HLut+fVrviW/8diC34wBRAHHKq/5B9gzy6z5l5k93lFZ928HnzjDC6RjJ1nzTwMAPfrwpJSGDexFSbFTVizNHmf+sfcca6ZIOuXoY5X2/OZj0okVKePvgAQAfVGQdvTRCmtGiVEfq0wrnbLn2G3Vv8opB6YSX/EjCqVPHcwvOgCi41MHpzVi8C+tj7Sy4s4ZZhOrdlOc7+jiI4y/NN+oC+fk8fGwACJlZJGjC+Yk+9h98ZF51r2s2aoBL0kzxqV0UkJP93x4UkrHTEzm9w4g2o6dmNKxCT1+nVSZ0oyx9n3v9u1I0rkz04k7VT+iUDp/FqfmAUTXF2cn71R9WbF07gw7j91WDvgknqrn1DyAqEviqXobT83nWDngpWSdqufUPIC4SNKpeltPzefYuzNJ581K6+AyO38z8kt5qaOL5th5egcABmLREWmVW/BpakE6uMzReTPtPnZbPeAL8xx96/g8lVvwqTxBmDBcumJunoYWxPP7A5BMQwscfWdunsYPN72TYJSXdM6mwjy7j91WD3jp/TvKvDxNiNkdZXSx9N15+TzuDiCWSoocXTkvX6Nj9oTpCcOl78yLRphZ82EzvWlqcfX9Jzu0rikS292nCSM6h3tZsf13EAAYjIYdrhY/2a6N20zvZPDKSzqDsyQiYRaZAS9JO9pc/fBvHVrZEJktd1FR2nnqinIHkBRbW1x976kOrdsS3WP3wWWdp+WjUO45kRrwktTS4erumoz+vDZreiv9NndySufPTkfqDgIAftjR5uoXyzJ66s3oHbtPqkzpvJlp6x9z31vkBnzO8vqsbl7aoYadpnfSu5Ii6aIj8nTkBOuf8gAAgXphY1a3Lu1QU4vpnfSurLjzde42vxRuXyI74CVpZ7ur/3zJ7pqfOzml82alNbwwWr/5AUBQtrW6urMmo79ZXPMnVaZ07oy0tW9i0xeRHvA5NtY81Q4A+/bChqxufdGumo96te8pFgNe6qz5/30jqz/VZowO+tKizt/8Pn4g1Q4AvdnW6uqxNzL6c21WWwwO+rLizs9zP+XAVKSrfU+xGfA5mayrZfWulqzOaPnb4X1rh45xNP+AtI6c4CgvFY87BwCEpSPr6oWNncfuVzaHd+yeOc7R/KlpzRrnKB2zY3fsBvye6re5+uOajP66Lqvtbf5ffnG+NG9KSvOnprX/yHjdMQDAlLqtncfuJ9dntbPd/8sfViB9pDylk6emNW54fI/dsR7wOe0ZV+ubXK3d4mpNo6u1ja7e2uoq04/vPOVIE0c4qhzV+b+KUkflpY4K0vG9cwCASa0dXY/ddc2usv04dqcdadJI77F7Somj/AQcuxMx4LvTlnH15vt3nKaWzl8CWjNSR1bKS0kFaakg7WhEoXbfIaL2GkgAiJs9h35za+exvG2PY3dhWspPOyop6jx2Ty5JbogldsADABBn0X8dAAAA6IIBDwBADDHgAQCIIQY8AAAxxIAHACCGGPAAAMQQAx4AgBhiwAMAEEMMeAAAYogBDwBADDHgAQCIIQY8AAAxxIAHACCGGPAAAMQQAx4AgBhiwAMAEEMMeAAAYogBDwBADDHgAQCIIQY8AAAxxIAHACCGGPAAAMQQAx4AgBhiwAMAEEMMeAAAYogBDwBADDHgAQCIIQY8AAAxxIAHACCGGPAAAMTQ/wGqSZWpymN+AAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 660x660 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(2.2, 2.2))\n",
    "G.plot(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "PZ = pd.DataFrame({'Z': [0, 1], 'P(Z)': ['1-r', 'r']})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Z</th>\n",
       "      <th>P(Z)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1-r</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>r</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Z P(Z)\n",
       "0  0  1-r\n",
       "1  1    r"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PZ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "MedicationProb = pd.DataFrame({\n",
    "    'X': [1, 1, 0, 0],\n",
    "    'Z': [0, 1, 0, 1],\n",
    "    ' P(X | Z)': [r'$q_{0}$', r'$q_{1}$', r'$1-q_0$', r'$1-q_{1}$']\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X</th>\n",
       "      <th>Z</th>\n",
       "      <th>P(X | Z)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>$q_{0}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>$q_{1}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>$1-q_0$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>$1-q_{1}$</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   X  Z   P(X | Z)\n",
       "0  1  0    $q_{0}$\n",
       "1  1  1    $q_{1}$\n",
       "2  0  0    $1-q_0$\n",
       "3  0  1  $1-q_{1}$"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MedicationProb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "DeathProb = pd.DataFrame({\n",
    "    'Y':[1, 1, 1, 1, 0, 0, 0, 0],\n",
    "    'X':[1, 0, 0, 1, 1, 0, 0, 1],\n",
    "    'Z':[0, 0, 1, 1, 0, 0, 1, 1],\n",
    "    'P(Y | X,Z)':[r'$p_{10}$', r'$p_{00}$', r'$p_{01}$', r'$p_{11}$', r'$1-p_{10}$', r'$1-p_{00}$', r'$1-p_{01}$', r'$1-p_{11}$']})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Y</th>\n",
       "      <th>X</th>\n",
       "      <th>Z</th>\n",
       "      <th>P(Y | X,Z)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>$p_{10}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>$p_{00}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>$p_{01}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>$p_{11}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>$1-p_{10}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>$1-p_{00}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>$1-p_{01}$</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>$1-p_{11}$</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Y  X  Z  P(Y | X,Z)\n",
       "0  1  1  0    $p_{10}$\n",
       "1  1  0  0    $p_{00}$\n",
       "2  1  0  1    $p_{01}$\n",
       "3  1  1  1    $p_{11}$\n",
       "4  0  1  0  $1-p_{10}$\n",
       "5  0  0  0  $1-p_{00}$\n",
       "6  0  0  1  $1-p_{01}$\n",
       "7  0  1  1  $1-p_{11}$"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DeathProb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"width: 100%; overflow: hidden;\">\n",
    "     <img src=\"data/D4Sci_logo_full.png\" alt=\"Data For Science, Inc\" align=\"center\" border=\"0\" width=300px> \n",
    "</div>"
   ]
  }
 ],
 "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.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}