Regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the OLS community, it is often believed that bivariate correlations are able to tell the existence of confounders. For instance, if Y $\\sim$ X, X $\\sim$ Z, and Y $\\sim$ Z, then controlling for Z would undermine the linear relationship between Y and X. Is it really so?  Using simulations, this notebook tries to debunk this myth and demonstrates the effect of control in OLS regressions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's draw a sample from a multivariate normal distribution $\\begin{bmatrix}\n",
    "X \\\\\n",
    "Y \\\\\n",
    "Z \\\\\n",
    "\\end{bmatrix} \\sim \\mathcal{N}(\\mathbf{\\mu}, \\mathbf{\\Sigma})$ with\n",
    "\n",
    "$\\mathbf{\\mu} = \\begin{bmatrix}\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "\\end{bmatrix}$ and $\\mathbf{\\Sigma} = \\begin{bmatrix}\n",
    "1 & i & j \\\\\n",
    "i & 1 & k \\\\\n",
    "j & k & 1 \\\\\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "1 & 0.9 & 0.9 \\\\\n",
    "0.9 & 1 & 0.9 \\\\\n",
    "0.9 & 0.9 & 1 \\\\\n",
    "\\end{bmatrix}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = np.array([0, 0, 0])\n",
    "cov = np.array([[1, 0.9, 0.9], [0.9, 1, 0.9], [0.9, 0.9, 1]])\n",
    "sample = np.random.multivariate_normal(mean, cov, 10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the covariance matrix has to be [positive semidefinite](https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Non-degenerate_case).\n",
    "We can verify whether a matrix is positive semidefinite by checking whether [all of its eigenvalues are non-negative](https://en.wikipedia.org/wiki/Positive-definite_matrix#Characterizations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.1 2.8 0.1]\n"
     ]
    }
   ],
   "source": [
    "eigenvalues = np.linalg.eig(cov)[0]\n",
    "print(eigenvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get our sample x, y, and z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sample[:,0]\n",
    "y = sample[:,1]\n",
    "z = sample[:,2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot y against z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABG+klEQVR4nO3df3QU9b038PdsSDY/SBaSCBskSAxUiRERlJKCVn6EJ0gR8epzrdVatd4HBI9I60W0itRrI23vxZ6iVK0Xep8UsVYBqTaPUBAVoYAhQkyrkAZDSaImgU0IZAO78/wRZtndzO7M7O7szO68X+fkHJPsj+/Cynz2+/38EERRFEFERERkAJvRCyAiIiLrYiBCREREhmEgQkRERIZhIEJERESGYSBCREREhmEgQkRERIZhIEJERESGYSBCREREhhlg9ALC8Xq9aG5uRnZ2NgRBMHo5REREpIIoiujq6sKwYcNgs4Xf8zB1INLc3IzCwkKjl0FEREQROHbsGIYPHx72NqYORLKzswH0vZCcnByDV0NERERqdHZ2orCw0HcdD8fUgYh0HJOTk8NAhIiIKMGoSatgsioREREZhoEIERERGYaBCBERERmGgQgREREZJm6BSGVlJQRBwOLFi+P1lERERGRycQlE9u3bh5deegljx46Nx9MRERFRgtA9EDl16hS+973v4eWXX8bgwYP1fjoiIiJKILoHIgsXLsTs2bMxY8YMxdu63W50dnYGfBEREVHy0rWh2YYNG1BTU4N9+/apun1lZSVWrFih55KIiIg083hF7G3swFddPRiSnY6JRblIsXEGWizoFogcO3YMDz30EN59912kp6erus+yZcuwZMkS3/dSi1giIiKjVNe1YMWWerS4enw/K3CkY/mcElSUFhi4suQgiKIo6vHAmzZtwrx585CSkuL7mcfjgSAIsNlscLvdAb+T09nZCYfDAZfLxRbvREQUd9V1LVhQVYPgC6W0F7LmzvEMRmRouX7rtiMyffp0HDp0KOBn99xzDy6//HIsXbpUMQghIiIykscrYsWW+n5BCACI6AtGVmypR3mJk8c0UdAtEMnOzkZpaWnAz7KyspCXl9fv50RERGazt7Ej4DgmmAigxdWDvY0dKCvOi9/Ckgw7qxIREcn4qit0EBLJ7UierlUzwd577714Ph0REVHEhmSrK7RQezuSxx0RIiIiGROLclHgSEeo7A8BfdUzE4ty47mspMNAhIiISEaKTcDyOSUA0C8Ykb5fPqeEiapRYiBCREQUQkVpAdbcOR5OR+Dxi9ORztLdGIlrjggREVGiqSgtQHmJk51VdcJAhIiISEGKTYhbia7V2skzECEion6sdjE0Cyu2k2cgQkREAax4MTSDUO3kW109WFBVk7Q5KUxWJSIiH+liGNxRVLoYVte1GLSy5KbUTh7oayfv8eoyHs5QDESIiAiAtS+GRtPSTj7ZMBAhIiIA1r4YGs3K7eQZiBAREQBrXwyNZuV28gxEiIgIgLUvhkazcjt5BiJERATA2hdDo1m5nTwDESIiAmDti6EZWLWdvCCKomnTnzs7O+FwOOByuZCTk2P0coiILIF9RIyVDM3ktFy/GYgQEVE/yXAxJONouX6zsyoREfUTz9kqZG3MESEiIiLDcEeEiIjIBKx6HMZAhIiIyGBWThDm0QwREZGBrD5okIEIERGRQThokIEIERGZkMcrYndDOzbXHsfuhvakvRBz0CBzRIiIyGSslC/BQYPcESEiIhNJ5nwJuV0eDhrkjggREZmEUr6EgL58ifISZ8KVtYba5XlidgkKHOlodfXIvm4BfbNmknnQIHdEiIjIFJI1XyLcLs/C9TW46aq+4yarDhpkIEJERKYQr3yJeCbCqqmKeeuTFjx/x9WWm7or4dEMERGZQjzyJeKdCKt2l2dwlh0fLp3GzqpERERGmXDJYORmpaKj+6zs76PNl5COSIJ3J6REWD12H7Ts8lh10CCPZoiIyHDVdS349i92hA1CgMjzJYxqHMaqGGUMRIiIyFChkjn9OTJSsXjGaJSXOCN6DqMSYScW5aLAkd4vEVUioO9oKJmrYpQwECEiIsOE26nwd/LMWazadhhTVm6PqJeIUY3DUmwCls8pAWDdqhglDESIiMgwSjsVwSJtbGbkEUlFaQHW3DneslUxSnRNVl2zZg3WrFmDo0ePAgCuuOIKPPnkk5g1a5aeT0tERAlC6w6EtHPy+MY6TLt8KNIGqPs8LR2RGNU4rKK0AOUlTktWxSjRdUdk+PDhePbZZ7F//37s378f06ZNw9y5c/Hpp5/q+bRERITEGBwX6Q5Ee3cvJlVuU70zYoYjEqkqZu64i1FWnMcg5DxBFMW4vjNzc3Pxi1/8Avfdd5/ibTs7O+FwOOByuZCTkxOH1RERJYdEGRzn8YqYsnJ7yJ0KJQKA5++4GoOz7Kp2GhLlzyXRabl+xy0Q8Xg8eP3113H33XfjwIEDKCkpUbwPAxEiIu1C9csQ0He0ce/kkSgvcZrmaEBaL4CIghGbAPhv9igFFh6vyCMSnZkqEDl06BDKysrQ09ODgQMHYv369bjxxhtlb+t2u+F2u33fd3Z2orCwkIEIEZFK0g6DmgRQM+0EyO1UREoKKZgIahwtgYjuVTOXXXYZamtrsWfPHixYsAB333036uvrZW9bWVkJh8Ph+yosLNR7eURESUVLFUqkFSh6qCgtwIdLp+HV+ydh0dTiqB5LPP+lR4Myij3dA5G0tDSMGjUK11xzDSorK3HVVVfhV7/6lextly1bBpfL5fs6duyY3ssjIkoqWqpQIu0oqlcSrJTMWXzRwJg8XiJO6rWiuM+aEUUx4PjFn91uh91uj/OKiIiSh9YqFP+OomrmnOid7Fld14Kn3/5b1I8jae2MbYMyij1dA5HHHnsMs2bNQmFhIbq6urBhwwa89957qK6u1vNpiYgsS6lfRihqdlL0HhoX6vGj0XFK/oMvmYeuRzNffvkl7rrrLlx22WWYPn06/vrXv6K6uhrl5eV6Pi0RkWWF65cRjtJOit5D49S2etcqNystxo9Isabrjsgrr7yi58MTESW9SEpNpZbiaqpQ1HYU1TI0Tssoe+n17TrytaaKmez0AejqOad4O6cjQ/VjkjHiniNCRETqRJOP4d9SfFt9K17ZdVT2diKAJ2aPUQxuYj00zuMVsXr7YazddRQnz5xVdR9/p1QEIVafapsoGIgQEZlQLPIxpCqUsuI8XFuUG3KH5Om3/wabTQj7eLEcGldd14JH3zyEk6e1ByASpSMcAZxqmyg4fZeIyGT0yMeoKC3AE7PlO1qr6SciJcGGuqwLULcDIQVY0QQhwXKzUgO+L+BUW1XMMouIOyJERCajRz6Gxyvi6bflm0mK6AskVmypR3mJM+Quwu3XjsCqbZ/3+7naoXF6JaTOHXcxZpY4I27ZbsWW72aaucNAhIjIZGKdjwFEF9wotV93hrmA+V/k27rcMWnhHuyt2mb85Pxuz97GDvzpYLPqgMJMF+R40bsMWysGIkREJhPLfAyJ2qBl15GvAy7gSr09Hp4xGoumjQ644EvBx7b6VmysPY6O7tgdw8hp7+7F6u2HsWHfMU0BhdkuyPGgdOynZmcs1pgjQkRkMrHKx/CnNmhZvaMBU1ZuR3Vdi+JRigBgw77AURzVdS2YsnI7vvvyHryy66juQYhk1bbD/XZbwuW+6N0Xxay07IzFCwMRIiKTUdOUTGtFyMSiXDhz1I3QaHH1YH5VDZb+8RNVF62n3qpD7zmvb4dBj+OXSEjD7x7f2Lc+f2a8IMeDHsd+0WIgQkRkQlJTMkdmar/fyf1MSYpNwHcnjtB0nz/WHFd1u/+7pwmXP/FnLPnDJzFPRI2F9u5eTKrcFrAzYsYLcjzocewXLQYiREQmJlfmevL0WcyvqsGvtn2u6ehgZH5WLJcWwCsCp3s9uj1+tDq6zwYc05jxghwPehz7RYuBCBGRCUk5DOGs2nYYk5/dHrb/hz+zXVTTBsT/EiTlfZjxghwP4Y791JZhxxoDESIiE1LKYZC0dio3I5NMuGSwKYbACehrQhact6E3/7wPM16Q40U69nM6AgNTp0GN4Fi+S0RkQlpzE5RKLqV+GR3dvbFYXlREAPPGXRxy/o3epD/bUMMBw/VFSRb+s4iMbuTGQISIKE60dPDUcoyi1GlVqRdIvA3OTMW0MUMNC0SOtp32/beZLsjxJs0iMhoDESKiONDawVPKYWh19agOIOR2UfRqqx6NE6fPAmLf6zei1HfDviYsmjbKF2yY5YJsVcwRISLSWaj+GlK/DrnqF/8cBrXkdlHU5prE2/+rb8Xt144ImSyqp2TsD5LIGIgQEelIzY5EqOoXX1KhQiOycBUeZu2D8T+7v8CqbZ/DkZmKQUF9UbLtKbo/v1n/XKyIRzNERDrSWv0SXLUg5TCs3n5E1eTb4DyU/IHquqkaxXX6LET0zawZmZ+FIdnpaO3swcOv1er6vGYrZbYyBiJERDqKtPoFQEBAsWjaKFzmHBi2wkMuD8WZY8egzFTZxmhmIA1a27DvGD5cOg0pNgHPbf0s4scT0Nd5Vgpw5H7vTML+IImMgQgRkY4iqX6RmySbm5WK/5hbig+XTpOt8AhVGfNlp9tUiapypNe9p6Ed+7/owHN/ORLR4wzKTMWzt1wJAFhQVQMBCHjtyd4fJFEJoiia9j3a2dkJh8MBl8uFnJwco5dDRKSZxytiysrtmqpfwpl62UWYMiofuQPtcOZc+GQ/ZeX2kEdAAoBMewq63eZtwQ4AjoxUuM5EvnPz+/u+icmj8wEoVylpKaUm7bRcvxmIEBHpTM8+HgWOdNx+bSFWbTuseNvs9AHo6jmnwyqMV+BI9x3tSEIFG1pLqUk7BiJERAaSuwBurW/FU299itZOt2HrKhycgWMnzhj2/Hr6jcrW5KGCQil8MaLFeTLScv1mjggRUQyF+7S969HpIatf4iEZg5DMtBT8n+sv9SX4hhOulFpKmlVqlU+xxz4iREQxEqpxWev5xmWrtx/GyPxMPDxjtGJvkFgSAGSk6t+bwwinez1Yte0wpqxUnkKsVErt3yqf4oc7IkRECtQkNip92gYQkMcxKCNV5pb6EAG4z5k7UTVarS75Piz+1JZSs9lZfDEQISIK452DLfjJ5rqAqbVyiY1aW6mfjKI6JBJe02YDxoaaoxW1pdRsdhZfPJohIgrhmbfr8cD6moAgBOjbvl9QVeM7CvB4Rew60mbEEsmP0tGKNEgwVPZHuFb5pB/uiBARyXjm7U/x8gdHQ/5eRN+nb69XxNNv/82Ug+WsKtTRijRIkM3OzIU7IkREQd452BI2CJG0uHrwwPoDDEJMJtzRim+QoCPwNk5HOkt3DcIdESKyFKXEU49XxE821xm4QoqU2jky0iBBdlY1BwYiRGQZatp+r9vV2C8nhMxP69FKik1AWXGevosiVRiIEJElhOqoKZV9/tv1RXjrkxYesyQoJ1u0JywGIkSU9NT0+Hjx/cZ4LolUsA+wwX3Oq3i7J2aPwQ8mF/FoJUExWZWIkpp03BLrnY7MVBucOaFLQSk6gzJTccfEEapum59tZxCSwHQNRCorK3HttdciOzsbQ4YMwc0334zPPvtMz6ckIvKprmvBlJXb8fTbf4v5Y58+68V3z18oeQmMvZOnz2L44AxVt2UDssSmayCyc+dOLFy4EHv27MHWrVtx7tw5zJw5E93d3Xo+LRFZmMcrYndDO3665VPMl5n7Ekvt3W4snjEaQ7LjNzcmGWSmqZt7k5uVFrYBGcAGZMlAEEUxbo1/v/76awwZMgQ7d+7E9ddfr3h7LWOEicjaPF4Rq7cfxtpdR+PePt0mJH8LdSO8ev8kuM70YkFVDQDI5vgMykzFs7dcySRVk9Fy/Y5rjojL5QIA5ObKR69utxudnZ0BX0RESqrrWjDhP7Zi1bbDcQ9CAAYhemk75fY1IHNkyg8JdJ0+G9BunxJP3AIRURSxZMkSTJkyBaWlpbK3qayshMPh8H0VFhbGa3lElKCkstyTp+MfgJC+nnrrU3i8IspLnEgfIH+5kmLAFVvq4WFEmJDiFogsWrQIBw8exKuvvhryNsuWLYPL5fJ9HTt2LF7LI6IEFK4slxJfe3cv9jZ2YG9jB1o73SFvpzTsjswtLn1EHnzwQbz11lt4//33MXz48JC3s9vtsNuZ9EVEF4Rryb63sYMNyJJcqAF20d6WzEPXQEQURTz44IPYuHEj3nvvPRQVFen5dESUZJRasm+rbzVwdRQPWkpzWcabmHQNRBYuXIj169dj8+bNyM7ORmtr3z8aDocDGRnq6sOJyJqUWrI/f8d4bKw9rvlxczNT0cF8koTgX5pb4EhHq6tH9hhO7bA7Middc0TWrFkDl8uFG264AQUFBb6v1157Tc+nJaIEp9SSXQTwk02H0NGtPaBgPqM5FOVlYvrlF4X8vYALA+xSbAKWzynx/Tz4doD6YXdkProfzRARaaUm9yPSXQ0jynupv8b202hsPw0AEATA/3JRIDPATirjDT6q47C7xMehd0RkuOCE1FbXGaOXRPF0Pgi5d/JIlJc4AxKS/VWUFqC8xBkyeZkSEwMRIjKUXEJqblaagSuieBPRd8Ty57pWPD77whGLf4Can2UHhL4mZ0Oy0/GdscMYgCQJBiJEZJhQCakd3b2GrIeM498LpKw4TzZA9Sd3fEOJiYEIEcVN8Cfcp95iMzIKtOtIG05092Lh+v4Bqj+pemrNneMZjCQ4BiJEpLu+gXRHsHZXI5NFKazVO47AJsgPuPMnHees2FKP8hInj2kSGAMRItJVdV0LHn3zEGfBkGpqS6yDj3MoMTEQISLdVNe1YP75Ee5EemFr98QWt6F3RGQtUlMyIr2xtXti444IESkKN3guFA6kI72xtXtyYCBCRGEpDZ4LhdvlpCe2dk8ePJohopCkPh/BOxtS6WR1XYvs/TxeEW1d7ngskUxuoD3yz7vhwgunI52lu0mCOyJEJEtp8Fyo0kmlRlRkLT//l7F4bFNkVVOOzFTc860ijMzP7NdZla3dkwcDESKSpZTjIVc6GapTKlmPdHwHIOLSbdfps3hu2+dYc+d4TB6dH8vlkYkwECEiWWpzPKTbhdtBoeQnALhl/MWYMvoiOHMuJJBO+I+tET8mm5ZZA3NEiEiW2pJI6XaskrE2EcAbNcfx8+q/w3WmFyk2Aau3H4m6kZ3/zhslJwYiRCRrYlEuChzpIRMGBfRtv0uffFklQ8CFROZ3DjZj7a7GmD0u31/Ji4EIEclKsQm+M/7gYESudJJNpQjo28EQASx982BM5wrx/ZW8GIgQUUgVpQVYc+d4OB2BFwG50klpB4UIALp6PKpu50gfAGeO+p03Sj5MViWisCpKC1Be4lTsrCrtoHC2DGlx75QiXObMxoKqGggInLrLpmXWwB0RIlKUYhNQVpyHueMuRllxXsiLQkVpAV6442oIvGaQCoMzU7Fo2mhNO2+UfLgjQkQxdePYYfjB0Q6s/egLo5dCMZSRasP/vmY4thxsQUd39LkfAoDKW670BbVqd94o+TAQIaKIhRqGN/OKAgYiSeZnt4yFMycdv9vdFPVjhZpVJO28kbUwECGyoEim6QYLNwyvvMSJAkc6+4okEWdOekxKaJ+YPQY/mFzEnQ7yYSBCZDGRTtMNfgy5Vu6trh7Mr6rBwzNGo+KKodwVSRJS1Uo0TcUE9OV8MAihYExWJbKQSKfpSjxeEbsOt+HRNw6FHIYHAKu2HWYQkkSkqhWlJnfhiGD1C8ljIEJkEUrTdIG+mR4er/y0mOq6FkxZuR3fe+WvMW1UReY1ODMVv/GrWgnX5E7JoMxUlJc4Y7xCSgYMRIgsQss03WChdlIoed08bhhWf3d8v+AhVKmtkpOnz3JeDMlijgiRRWidpivhVF1r2lTbjE21zQEJyP4JzjsfmYqPvziBr7p6cPjLLqze0aD4mJwXQ3IYiBBZhNZpuhJO1bU2KQF5UGZqwCRdKUCZO+5i7G5oVxWIcF4MyeHRDJFFqEk0HGhPQWtnD3Y3tPtyRfgp1tqknTD/IAQITHDWOqmZyB8DESKLkBINwx2xnHJ78PBrtfjuy3swZeV2VNe18FMsyfJPcAagaVIzkT8GIkQWUl7ixKDMVFW3lT7xnuju5VRdkuWf4Mx5MRQp5ogQWcjexo5+W+yhiOj7NPv02/X4zlgnXv7gqJ5LowQmHd9xXgxFgoEIkYVozfeQPvG+UXNcnwVRUvA/vuO8GNJK16OZ999/H3PmzMGwYcMgCAI2bdqk59MRkYJI8z1iMW2Vkg+TUCkWdA1Euru7cdVVV2H16tV6Pg0RqSRVNxBpxSRU0ouuRzOzZs3CrFmz9HwKItJAqpyZX1Wj+j4D7QNwyn1Ox1VRIhiclYaO7l7f906NgxKJQjFVjojb7Ybb7fZ939nZaeBqiJJTRWkBHp4xGqu2HVZ1e1FkT1UCnpg9Bk5HBpNQKeZMVb5bWVkJh8Ph+yosLDR6SURJadG00chMS1F12+5ej86roUTgdGSgrDgPc8ddjLLiPAYhFDOmCkSWLVsGl8vl+zp27JjRSyJKSlvrW3GaAQapwIRU0pupjmbsdjvsdrvRyyBKatIQOyIlTEileDDVjggR6Y9D7CiU3KzArrvsikrxoOuOyKlTp3DkyBHf942NjaitrUVubi5GjBih51MTWZrHK4bsbskhdiQnLysNu5dNx8dfnGBCKsWVroHI/v37MXXqVN/3S5YsAQDcfffdWLdunZ5PTWQ5UvCxtb4Vm2qbA0otC/xKLf/xdbeBqySzmjtuGNIG2NgVleJO10DkhhtuYOkfURTC7Wz432b19sNYu+soTp6R74Da4urB/KoaLJ4+Cr/76GgcVk6JprzEafQSyKJMlaxKRBdU17VgxZb6gHyOgqAmUtV1LXj0zUOqB9k995cjyjciy2FVDBmJgQiRCVXXtWBBVQ2C9xNbXT1YUFWDNXeOBwDZ2xCppaYqxn9XLj/LDghA2yk3c0goZhiIEJmMVF4rF2CI6Lt4rNhSD1EUGYSQJgIQ8J5RatMutyvnL3iHjigSDESITEapvFYEWH5LEREB3PnNEbDZBFySm4m7ykYibYB8F4dQu3L+/HfoGIxQpBiIEJkMy2tJT1V/bfL9928/bJTd0Qi3K+fPf4euvMTJYxqKCBuaEZnMkOx0o5dAFiHtaFTXtQT8XEvTO2mHbm9jhw4rJCvgjgiRyUwsykWBIx2trp6Qn0gHZaTizFkP3Oe8cV0bJRfp/fXYxkM4c9YLZ05fAmoku3LcyaNIMRAhMpkUm4Dlc0qwoKqmX3KhJFS/EKJIdHSfxcOv1QLoS0C9/Vrtna+5k0eR4tEMkQlVlBZgzZ3j4XTwH3eKr1ZXD57b9jkGZaZCTcYHp/NStBiIEJlURWkBdj4yFY/fOAaZqSlGL4cSyODMVOUbhSDtwEm7ceGCEU7npVjg0QyRSSn1cCAK5cnvlMDpyPA1IfvR65/gy87QOUfBRAAnTp/FwzNGY8O+YyHfg0p9SIjUYCBCZEJqejgQheJ0ZAQMr3vqpvA5R6GMzM/Ch0unsbMq6YqBCJHJqO3hQCRHLl9DyjnSusM2JDsdKTaBE3lJVwxEiExGSw8HsoYCRzrO9HpUVUvdfm2h7M8rSgtQXuLE3sYOtLrO4Om3/4YT3b2yAa+AvmMXJqBSPDBZlchk2I+Bgv3y1qvw7L9cCQHhk0cBYNW2w5iycnu/JmUAfLsb88YPx8/mlQIyj8cEVIo3BiJEJsN+DBSsrdutqaQ7VMdUf6Eez+lI5+wYiisezRDFmf9Y9eCEv95zXhw6dhIpAuBhkgidJwWncscrHd29/W6vdgaM/+PJvR+J4oGBCFEcyZXkSqPUDzSdwEsfNEJkAEJ+crNSA3I1pOOV3Q3tskGIxH8GTLhkUyajktEYiBDFSaiS3FZXD+ZX1RiyJjK/uyaNxJ8ONvfbrVCbS8ScIzI7BiJEcRCuJJcbIBRKZloKfvWXw77vC/waiKnNJWLOEZkdk1WJ4oAluRSJ072egO/9k1ClKc2hsjk4A4YSBQMRojjg9jjFgrR7tmJLPYC+EluAJbiU2BiIEMUBt8cpVvyTUEOV4OZmpeH5O1iCS4mBOSJEOvEv083PsmNQZipOnlbujEmkhrTLVlFaAK9XxE8216Gju+/91d7di6ffrofNBgYjZHoMRMgSwvXu0AMn55LepF226roWLFx/QLYaa0FVDZuTkekxEKGkF653hx7/QHNyLunJfw6MUjWWmqZmREZjjgglNSkoCN6ZUNMCOxKcnEt6Ck5CVarG8s8nITIr7ohQ0orVp0W5Yx0Askc9LNOlWMqyp6DbfaGE1xm0k8emZpQMGIhQ0tLyaTFUi2u5Y51BmakAEJB4Kh31uM95Y7N4IgDdbg9ys9Jw87hhKC9x9sttYlMzSgY8mqGkFe2nxVDHOidPn+1X/SId9RxtOx3ZYolCONHdi7W7jsJ1prffzh2bmlEyYCBCSSuaT4tacz2k223Y1wRnjj3khYFIK/8mZh5v4DsyxSawqRklPAYilFA8XhG7G9qxufY4dje09/uH2V80nxYjyfWQjnquGu5gsirFVLik01BNzZyOdJbuUkJgjgglDK1luNKnxQVVNRAQOFxO6dNiNMl9/6/+q4jvS8kpMy2l39yYSIR6X1aUFqC8xBnXXjlEscIdEUoIkZbhRvppkcl9FEs3ljpj8jjh3pcpNgFlxXmYO+5ilBXnMQihhBGXHZEXXngBv/jFL9DS0oIrrrgCzz33HK677rp4PDUlgWjLcCP5tCgd67S6enjMQmHZBwhwnwv9LhmUmYqf3TIWuxraw76fbAIgipD9vX8TM6Jko/uOyGuvvYbFixfj8ccfx4EDB3Dddddh1qxZaGpq0vupKUnEommT1k+L4ZIAiQIJEBD6ffLsLVcibYAtbFKpAOD+64pC/h5g0iklL90Dkf/6r//Cfffdhx/+8IcYM2YMnnvuORQWFmLNmjV6PzUliXg1bQpOhC0vccoe6xD5c5/zYvGMb/R7nxQ40vEbv+M/pWPCZTeWMOmULEnXo5ne3l58/PHHePTRRwN+PnPmTHz00Uf9bu92u+F2u33fd3Z26rk8ShBq8zXys+zY3dAeUbJeuETYD5dOw55/tOOB39fAdYbTc6m/sx4P/v1/XYaO7l7kDrTDmSP//lM6JmTSKVmRroFIW1sbPB4Phg4dGvDzoUOHorW1td/tKysrsWLFCj2XRAlIKV9DAODITMWPXv8ErZ3aB9uFGlLX4je9tA+zRUje6h0Nvv+W3nehggfpmDAUpd8TJZu4VM0IQuD/kKIo9vsZACxbtgwul8v3dezYsXgsj0xOqWmTiL5up/5BCKBusJ1S4zIRwI9e/wTzq2rgOnMu0pdAFqLXQEWiZKVrIJKfn4+UlJR+ux9fffVVv10SALDb7cjJyQn4IgJCn68PzbH7Zr8EC9eRUqKmcZn/0DEiJWred0R0ga5HM2lpaZgwYQK2bt2KefPm+X6+detWzJ07V8+npiQkd37uFUV877d/DXkfpcF2nEpKelAzUJGI+ujeR2TJkiW46667cM0116CsrAwvvfQSmpqaMH/+fL2fmpJQ8Pn55trjqu4XHHB4vCL2Nnbg8JddMV0fWcNA+wCccisf1THQJVKmeyDyr//6r2hvb8dPf/pTtLS0oLS0FO+88w4uueQSvZ+aLCCSwXZyFTJEWqgJQgB26CVSIy6dVR944AE88MAD8Xgqshg1FTX+HSlDVcgQaTUoIxWuM2fZCZUoSpw1QwlNyxh0pQoZIi3umTwSADuhEkWLgQglPDWD7TxeEet2NfI4hqImoK9XyKJpo9kJlSgG4nI0Q6S3cB0pmRNCsSQCuOmqAqTYBHZCJYoBBiKUFKQqmOCLwTsHm/HA+gNGL4+SzEvvN+LqEYNRUVrATqhEUWIgQgkv1JyY74wtwCsfNhq4MkpmK7bUo7zEyd0PoigxR4QSmlQFE3zs0uLqwcsfNIKNLUlObpZ8N161/BuWEVF0uCNCphbqyEX6HatgKBKpNgGDMlNx8nR005TZsIwoegxEyLRCHblIE3XVzIkhkvNVV29MAlg2LCOKHo9myJRCHbn4Tzblp1GKVLRBiFTCy4ZlRNFjIEKm4vGK2HWkDY++cUj2YuE/2TR/oD2eSyMCwIZlRLHGoxkyDbX9PqREQYgI296dSA9Ov+NBIooeAxEylJSMurW+Ff+966im+7Z1u7F8TgkWVNVAQPTb7UThDMpIxfPfG49Jl+ZxJ4QohhiIkGGi7Xg6JDsdZcV5WHPneDz65qGoKyCI5Eghx7P/ciUmj8o3dC1EyYiBCBkimim40mTTCZcMxu6Gduxt7GAQQhFT2k3jUQyRvhiIUNxF0/9D+nR601UF+PYvdrB8l6Ly8IxvYMO+pn4l4rdfOwIj8zM5O4YoDhiIUNxF0/9jaI4dc8cNw0vvNzInhKLSN0F3FBZNG8WhdUQGYiBCcRdd/w8Br+3/J4MQiphc+S2H1hEZh4EIxV003ShbO3kUQ9FhzgeRuTAQobibcMlg2ARwIB3F3eM3jsG9U4p49EJkIuysSnH38RcnGISQIUoKchiEEJkMAxGKO86IIaO0dbuNXgIRBWEgQnHHiaVkFL73iMyHgQjF3cSiXBQ40sENcooXTsslMi8GIhR3KTYBy+eUGL0MsghOyyUyNwYiZIiK0gIsnvENo5dBSUQAMDgzFc4ce8DPnY50rLlzPMt1iUyK5btkmJH5mUYvgRKEUrm3tM9RecuVKC9xslMqUQJhIEKGOdp22uglUALIzUrFrqXTUXvsJL7q6sHRttN4dW9TQHO74CZl7JRKlDgYiJAhPF4Rr+5tMnoZZHICgJ/NuxIZaSkBwQXnwxAlDwYiFFcer4i9jR3YdeRrtmunsArCtGJPsQnc9SBKEgxESDdS0CF9aj3R7cZP/1SP1k42laJA2ekpuHV8IYYPzkBuVhqcjgzuchBZBAMR0kV1XQtWbKlHi4u7HqSsq8eDdR8dZXULkQWxfJdirrquBQuqahiEkCYigBVb6uHhICIiS2EgQjHl8YpYsaUevJRQJFpcPdjb2GH0MogojhiIUEztbezgTghFZWt9a8D3Hq+I3Q3t2Fx7HLsb2rljQpRkdM0ReeaZZ/D222+jtrYWaWlpOHnypJ5PRzEQnGCqNWGQlTAUrc21zXh8dl87drlco3DVNESUeHQNRHp7e3HbbbehrKwMr7zyip5PRTEQi3/0O06xIoai097di72NHXCd6cWCqpp+x3ytrh4sqKphYitRktD1aGbFihV4+OGHceWVV+r5NBQDoRJMpX/0q+taVD1OblaaHssji2nt7AmZayT9jImtRMnBVDkibrcbnZ2dAV+kv3AJplr/0Xc6MmK6NrKmjlPusLlGIpjYSpQsTBWIVFZWwuFw+L4KCwuNXpIlKCWYavlHf2JRLgoc6TFcHVmJgL7jQLU7a191MSeJKNFpDkSeeuopCIIQ9mv//v0RLWbZsmVwuVy+r2PHjkX0OKSexyti15GvVd1WzT/6KTYBy+eUgP0wSUnwe0T6fvmcEtU7a0OyGfQSJTrNyaqLFi3C7bffHvY2I0eOjGgxdrsddrs9ovuSelJlzNb6VmyqbUZHd6+q+6n9R7+8xInFM0Zj7a6jOHnmbDRLpSR16/iLsauhPWAnzn+CrscrosCRjlZXj+yRoXD+9hOLcuO2ZiLSh+ZAJD8/H/n5+XqsheIgktbrWv7RZ2t3UmITgJ/dMhYpNiFkqbi0s7agqgYCEBCM+O+ccBYNUeLTtXy3qakJHR0daGpqgsfjQW1tLQBg1KhRGDhwoJ5PTTKkyhgtdQZa/tGP5PHJeqaPGYK0AX2nwuEm6FaUFmDNneP7BbZO9hEhSiq6BiJPPvkkfve73/m+v/rqqwEAO3bswA033KDnU1OQSFuvq/1Hn63dSa26453weEVVuxkVpQUoL3FG1WSPiMxN10Bk3bp1WLdunZ5PQSpF0nr9gRuKMXlUPtpOubHrSBsgAm3dbtmLAVu7k1pSBVa43RB/KTZB9W2JKPHoGoiQeURS5rhh3zG88F6D7O+y0wfg1gnDMbPEiYlFudgWNB+EKByW3RKRhIGIRURS5hiumqar5xzW7jqKtbuOIistBd29nmiWRxbDslsikpiqoRnpR2o0psfJOoMQkhQ40uHMsYd8n0kNy1h2S0QSBiIWIZVDEulp+ZwSPHXTFQDCNyxjsikRSRiIWIhUDjkoI9XopZDJRPuesAnAC3f0TcOV3mfOoFb/Tkc6J+YSUT/MEbGYitICZKen4nu//avRSyGTEADcM3kkVm07HPFjrP7u1bhx7IUAg2W3RKQWAxELmnRpXtj22WQdzhw7nrrpCpSXOLFh3zHFEuyMVBvOnPX6vi8I02eGZbdEpAYDEQsK1z6brOPhGd/AommjfLsU0nsi1Pvh/1xfhH+vGBN2l0OaY8RdECJSSxBF0bTXoc7OTjgcDrhcLuTk5Bi9nKTzzsEW/GRzneqhd5Q8bh0/HCtvHdsviFi9/QjW7moMGFaYl5WGp+eWBhy9yJGbMxRux4SIkpeW6zcDEYuSu2hk2VPQ7WYprlX4Bwly74dBGam4Z/JILJo2OuI5Q9K9mKRKZC1art+smrEg6aIRnA8gBSFC0DVHADDhkkHxWRzFTaurBwuqalD5Tr3s+8F15iye23YYWxW65oabMyT9bMWWeni8pv3MQ0QGYiBiER6viN0N7dh44Dge23gobF6I3B7Zx1+cRFZaim7ro/gTz3+9/EFjVEGE0pwhERfmyxARBWOyqgXIbbtrIV2C2EE1OYXbqPAPIkJVwKidG8P5MkQkh4FIkgt1dk+kRbggQu3cGM6XISI5PJpJYuHO7om0CBdEKM0x4nwZIgqHgUgS2/OP9oiPY8g6bEL/uTASNUGE/xwjzpchIq0YiCQhj1fEr7Z9jvt/t9/opZCJCee/7r+uyPd98O8BdUEE58sQUaSYI5Jkquta8Oibh3Dy9FnlG6vE7qvJyenXR+TqEYP7JTQ7NTYj43wZIooEG5olkWgTU6WAY1BmakAg48gYANeZc7FYIpnAoMxUPP/d8ZhUnMf27ESkCy3Xb+6IJIlYJKYOzbHjuxNHYEReFjpOuZGblYYh2elY9GpNzNZJxpFCimdvuRKTR+f3+z2H1BGRERiIJAmlplJygiepftXlDhgFX+BIx+3XjsCJGB7zkHG0HrUQEcUDA5EE5r+VfvjLLs339w9CgP6NrVpcPVi17fNolkg6+tXt47BiS73i0MKpl12Ef7u+mEctRGRKDEQSRPD5/YluN55++28sz7WgwZmpqLzlSlSUFsA+wBY2L6i8ZAhe/v61cV0fEZEWDEQSQLQt2il5LJxajCXll/l2NqSy2eD3R076APzs5ivxnXHDjFoqEZEqDERMji3ayd+UURf1O15h2SwRJTIGIibGFu0kEQDkZqWh1XUGuxva+wUarHghokTFQMTEIqmEkdw8bhg21TbHeEVkFBFAe3cvHv7DJwD6KppYAUNEyYAt3g3i8YrY3dCOzbXHsbuhHR6ZWezRjE0fPjgzmuWRybW6erCgqgbVdS1GL4WIKCrcETGAXPKp3Cfco22nNT+2gL5+EWXFeVi940gslksGCtVeXzz/uxVb6lFe4mQ+CBElLO6IxJmUfBp85CJ9wn3nYDN2N7Tj6S2fRtTDQwTwxOwxmHRpHgZlpGq+/9WFDs33If2Eyw8Scb7Xy9bPQ+6qERGZHWfNxJHHK2LKyu1h8z5sQv/GYlpJuyuftZ7SHMzkZaWhXaFBFpkT80aIyCy0XL+5IxJHapJPtQQhs690yv5c2l3pOqO9NXt7dy/SU/m2SETMGyGiRMQrThxFk3wq58Mj7bI/F89/vbKrMaLH7Qlq/U7GEYQLw+qUSDHsii31PKYhooTBQCSOhmSnx/TxXAo7HrwUJT5RvJCYqur26Msb2dvYoeOqiIhiR7dA5OjRo7jvvvtQVFSEjIwMFBcXY/ny5ejttW7+wcSiXBQ40lVfVIgEAIMyUzE0R1sQG+vdNyIivehWvvv3v/8dXq8XL774IkaNGoW6ujrcf//96O7uxi9/+Uu9ntbUUmwCls8pwYKqmpBlmUT+RAAnT5/F7+8bD5tNwK4jX2P1jgbF+8V6942ISC+67YhUVFRg7dq1mDlzJi699FLcdNNN+PGPf4w333xTr6dMCNKQMqcj8EKhtg0E20WYy+Lpo5GblRbws0GZqRiUqb10Opy2bjfKivPwcPllYXfVBPRVz0wsyo3p8xMR6SWuDc1cLhdyc/kPpNyQshPdvXhgfY3ifZmDaC5FF2Vh3+Mz+g2cA/qqpNTuYCiRdjjC7apJwcnyOSVscEZECSNugUhDQwN+/etf4z//8z9D3sbtdsPtdvu+7+zsjMfSDCE3pOzeL0biv3cdNWZBFJH8LHvIgXNlxXmYWJSLN2qOo9XVI3sUJ6CvMiZUgCl1yvXf4ZB21YK78zrZR4SIEpDmo5mnnnoKgiCE/dq/f3/AfZqbm1FRUYHbbrsNP/zhD0M+dmVlJRwOh++rsLBQ+ytKYOUl8n1BguVmpUaU8MrPyLH3o9c/Cdu3Q9rBAPr/+Uvf339dUV9AEuL3cjscFaUF+HDpNLx6/yT86vZxePX+Sfhw6TQGIUSUcDR3Vm1ra0NbW1vY24wcORLp6X1byc3NzZg6dSq++c1vYt26dbDZQsc+cjsihYWFSdNZVYnUeTXcp2enIx1PzB6DhesPAGDCq9Gk8GDNnePDBgFK84XUzh8iIkoEWjqr6tri/fjx45g6dSomTJiAqqoqpKSkaLp/srV4V0OaRQPIn/9LFzy5CxcZQwoQP1w6LWxuhscrYk9DO3b/ow1A33HOpEvzfPfxeMV+uSbM9SCiRGSKQKS5uRnf/va3MWLECPzP//xPQBDidKo7grBiIAIof3qWLlitnT3Ydfhr/LHmuIGrJcmr90+SzRWRcNeDiKxCy/Vbt2TVd999F0eOHMGRI0cwfPjwgN+ZeM6eKchV1UifjrkTYl7hmohJO13B73xpPozS0Q4RUbLSrY/ID37wA4iiKPtFyqRKjLnjLkZZcZ4vCFlQVcMgRKXUOB9rhGoi5vGKWLGlXjafh/NhiMjqLD9rxuMVsbuhHZtrj2N3Q3vEF4NYPU64xw91MSN5Z70i5ox1wpGub5W6UhMxpanLnA9DRFYW14ZmZhOrM/t4nP0rXcxI3l/+/jXsA/SLt9U0EVM794XzYYjIiiy7IxLqmEM6sw/XGyLax4lk94QXqcic7vXgxOnwU4rlhOrpEdy63elIV8zvUDv3hfNhiMiKLLkjonRmL6DvzL68xKlYjqn1cSLdPTH6InX96Hy8fzh8/5hkkZVmAyCgu9fj+5nUtTRUEnE40tRlpf4wnA9DRFZkyR2RWJ3Za32caHZhpIuZUV0lZo+1TkVHd6/XF4QMykjFwzNG+7qWyiURK1HTXZXzYYjIqiwZiMTqzF7L40RbOeF/MTPCtvovDXtuI7nOnMVz2w5ja31rVI8TauqymqMdIqJkZsmjmVid2Wt5HLW7J3v+0Q6bIMhu/UsXs0ffOISTZ7TnPUTDqomyWo7qlITrD0NEZFWWDERidWav5XH+dLBZ1doW/r4mIMgYlJGKeyYXYdG0UUixCagoLUBW2gDc9d97VT1erNQ1J+8kZCX+R2zhOqeqEWpSLxGRVVnyaCZWZ/bhHkciPY7a3ZPgnY6TZ85i1bbPMeE/tqK6rgXVdS145I+fqHosii1WLhERxZ4lAxEgdmf20uM4gso6AQT8LNpk05Onz2J+VQ3mV9WgtdOtfAcKqcCRjhfuGK/57yM/y67bmoiIrMqSRzOSWJ7Zn5TpVeE6fTZgjsjyOSVYUFUDAdCtQ6qej50McrNSsfORqUgbYIPNBt+kY1WYykFEFHOW3RGRRFKO6U+qhpETXA0TahcmuElWNBiEhNfRfRYff3ECwIXdrEEZ6v78205xJ4qIKNYsvSMSC1p6iZQV58nuwni9Ir73yl/jt2iLa+288PdVUVqA7PRUfO+3yn/+RjeVIyJKRgxEohRJT5Lgyonec14eqcTR03/6FBmpNl8e0KRL89j5lIjIIJY/molWLHqSrHmvgUFIHHV0nw3oZMvOp0RExmEgokK4IXVK1TBKI+Kr61qwatvnsV+0xakJGfw72bLzKRGRMXg0o0BpSJ30aVquGkbp03S4RFfSTjpCeWL2GPxk86fo6O4NeVu5JmXsfEpEFH/cEQlD7ZC6SD9NKyW6knbL55TgxrHD8MTsMapuH5zjE20VFRERacMdkRCUhtQFzx+J5NN0JJ06pUf7t+uL8PrHx8N+6g+WlZYSMNpeq/KSIag73tlvd+jxWWPw6MZDOOU+F/FjRysvKw3PzCv1BX1OR4aq+7EShojIWAxEZHi8ItbtatRUlgtonyMSyUXQ6Xcs9KOZl2NS5TZ0dKsbgJedPgAvff8a/HH/MWysVTf7xl/d8U7sfGQqPv7iRECwtbW+1dAgJMuegt3LpiNtwIUNvljNEyIiIn0xEAkilxMSTjTzR5QulgDgzLHjF7dehb82dgAQUXZpPiadD3bSBtjws3lXYr7K7qCtnW78vaUTWemR/bW3uHrw8RcnAoItM+S5pKXY+u08RZO7Q0RE8cNAxI+UE6KllFbNrobHK8oe2ai5WM4dNwz//sZBX2C0ekdDQLIs0NeZVa7FvJyn3/6bhlfXX3DgFes8l0j6qZw4fVZ2Mq6UuxMcWDqD/vyIiMg4DETOC5cTIkft1r5S1U24i+VNVxXgpfcb+62pxdWD+VU1uP+6Ivz2g/6/11Nw4LWtvjWmjz84K01T3osk1M4UK2GIiMyNgch5Wj7Zq93aD7XDIlXdSBU1chfLcYWDMHnlX8IGGS9/0KhqvbGSkz4gIPDyeEVsrD0e0+d4YvYYNHWcxqpthzXdL9zOlNbcHSIiih+W756nJddDTZMrpaobEcBjGw+h95wXQGDZqOtMLyav3K46CTVerh4xOCDw2tvYoWqNA+0DkJWWouo5hmSnY8O+Y6rXpNQwjoiIzI2ByHlqK1iemD0GHy6dpphfoGaHpaP7LCZV/sXXjwS4sIsSyfGE3q4fnR/wvdrg7brR+Xh6binSU0O/3aSAAgI05ZyI6Ps74VELEVFiYiByntpW7T+YXKTqoqf2It3R3etrjqY1TyWeBADfGJod0N5ebfD257pWLHn9E/Sc9YZ8bKDvqKvtlFvz2p5++28BwRwRESUOBiLnxXrwmZYeISKAxzfW4aMjbabttCoCuOu/92LKyu2+i75S8KaW/1FXJL1VgjvdEhFR4mAg4ieWg88mFuViUGaq6tu3d/diftXHqm9vFP+LfrjgLRwBfZ1QV/3vq/Dq/ZMCjromXDIYuVlpmtYk7dH4D7EjIqLEwKqZIFrKPUP1B4lUNO3X/RU4+qpu/lwX29JaoH97+1Dlx0qP0d7dC6cjI6CaRSp1jiQ/Rq7TLRERmR8DERlqyj2V+oPsbexQ3WQsVu6bPBIzSpyYWJSLvY0dugQiQP+LfnDwdvjLLqze0aD4OP55NO8cbMYD6w9EvbZoOt0SEVH88WgmAmqm8sb7gigAeKeu1bcrE6v8jXD8X6N/+fHkURepur+UD/LOwRYsejV8EBKu4kbuMYmIKDEwENFIqT8I0HdskZ9lj+p5BmUE5pfkKeRN+O9SAIHJt3oJddFXW4E0sSgX1XUteGB9DZRSO0JV3Mg9JhERJQ4ezWik1B9ECgggQHGgXTjP3zEeNpvgyz9pdZ3Bw3/4RPF+/rsUUv7GYxvrYtqXRKm9vdqBcwBiMjCPQ+yIiBIXd0Q0Unvk0nbKHXFFSYEjHZOK83xHHWXFeXA6MlTdP3iXoqK0AHuWTUdulvoKHqX1AcoXfTUVSLEamBdJVRMREZmDrjsiN910E2pra/HVV19h8ODBmDFjBlauXIlhw4bp+bS6UpuDMCQ7HWXFebIVJYMzU3Hi9FlN4+ml445QOyzhdinSBtjws3lXYkFVDQDt0239OTJTcc+3ilBe4lS8rVIFUrR5NIumjsLkUfkcYkdElMB0DUSmTp2Kxx57DAUFBTh+/Dh+/OMf49Zbb8VHH32k59PqQirVbXWdQW5WGk5096oKCEJdjLfWt2oaT6/2uCPUBbmitAD/dn0RXv6gEaLfnW0CcOOVBSgvGYoh2ek40e3GT/9Uj9bOCx1Os9JsECHgdK8HJ0+fxaptn2PDvqaQa/X/85Je83fGDuu3tmgTS0cPHchSXSKiBCeIohi3DlBvvfUWbr75ZrjdbqSmKh8VdHZ2wuFwwOVyIScnJw4rlCdXqitHusyqPSaIpA+JUtlwuPvJTQKW1u2/Zv91HW07jee2fd7vfuFeq9o1erwipqzcHnEezav3T2IgQkRkQlqu33ELRDo6OrBgwQIcP34cH374oext3G433O4Ln8Q7OztRWFhoaCAS7gIeTE1AEAtaAxjpgh8qkJJ2cT5cOi3gcSK5X6g/r1CBi3R7QP2RUaj1EhGROWgJRHRPVl26dCmysrKQl5eHpqYmbN68OeRtKysr4XA4fF+FhYV6L8/H4xWxu6Edm2uPY3dDOzxeUdUQutysVKz613H9WpXryb9nR1lxnuLFWG2lj1T6G+n91JY2+7dhD5XUOvh8e/xYzP0hIiLz0pwj8tRTT2HFihVhb7Nv3z5cc801AIBHHnkE9913H7744gusWLEC3//+9/GnP/0JgtD/IrJs2TIsWbLE9720I6K3UEcJt19bqHgc09F9Fs6cdFMfEahNCg2+ndb7aQlc/P+8YpVHQ0REiUdzILJo0SLcfvvtYW8zcuRI33/n5+cjPz8f3/jGNzBmzBgUFhZiz549KCsr63c/u90Ouz26RmBahTpKaHX1YNW2w6oew+xtxbVU+kRzv0gDHkC+rb6WuT9ERJSYNAciUmARCSkdxT8PxEhqjhLUMHtb8UhLf7XeL9KAJxw1c3+IiChx6ZYjsnfvXqxevRq1tbX44osvsGPHDtxxxx0oLi6W3Q0xQrQNtRKlrbh/u3ctORda76eltTsRERGgYyCSkZGBN998E9OnT8dll12Ge++9F6Wlpdi5c2fcj19C0XKkkuhJk2o6nUZ7v0gDHiIisq649hHRSu8+Irsb2vHdl/co3u7hGd/Ahn1Nmnt3mFEkvUu03i/SXidERJQcTNlHJBJ6ByJKDbX8+1UASNqkyUiDk3g/JhERJQYt129LT9/V2jY9GZMm9dq9YJIpERGpYfnpu5HmTiQDqXQ5OGG31dWDBVU1qK5rMWhlRERkFZbeEZFYsV+FUumygL4uqOUlzqT+cyAiImMxEDnPakcJkXZBJSIiiiXLH81YVTRdUImIiGKFgYhF6dEFlYiISCsGIhbFLqhERGQGlgxEPF4Ruxvasbn2OHY3tAeMpbcKdkElIiIzsFyyKrt+XiCVLgf/eTgt+udBRETxZ6nOqlLfjOAXLH3mT/a+IaGwCyoREcUSO6vKYN+M0KxWukxEROZhmRwRLX0zzIh5LURElIwssyOSyH0zmNdCRETJyjI7IonaN4PzYIiIKJlZJhBJxL4ZSnktQF9eC49piIgoUVkmEEnEvhmJntdCRESkxDKBCHChb4bTEXj84nSkm7J0N5HzWoiIiNSwTLKqpKK0AOUlzoTom5GoeS1ERERqWS4QARKnb4aU19Lq6pHNExHQt5tjprwWIiIiLSx1NJNoEjGvhYiISAsGIiaXaHktREREWljyaCbRJFJeCxERkRYMRBJEouS1EBERacGjGSIiIjIMAxEiIiIyDAMRIiIiMgwDESIiIjIMAxEiIiIyDAMRIiIiMgwDESIiIjIMAxEiIiIyDAMRIiIiMoypO6uKYt/M2c7OToNXQkRERGpJ123pOh6OqQORrq4uAEBhYaHBKyEiIiKturq64HA4wt5GENWEKwbxer1obm5GdnY2BCE5Brx1dnaisLAQx44dQ05OjtHL0RVfa/KxyusE+FqTkVVeJ2D8axVFEV1dXRg2bBhstvBZIKbeEbHZbBg+fLjRy9BFTk5O0v+PIOFrTT5WeZ0AX2syssrrBIx9rUo7IRImqxIREZFhGIgQERGRYRiIxJndbsfy5ctht9uNXoru+FqTj1VeJ8DXmoys8jqBxHqtpk5WJSIiouTGHREiIiIyDAMRIiIiMgwDESIiIjIMAxEiIiIyDAMRg910000YMWIE0tPTUVBQgLvuugvNzc1GLyumjh49ivvuuw9FRUXIyMhAcXExli9fjt7eXqOXpotnnnkG3/rWt5CZmYlBgwYZvZyYeuGFF1BUVIT09HRMmDABH3zwgdFLirn3338fc+bMwbBhwyAIAjZt2mT0knRRWVmJa6+9FtnZ2RgyZAhuvvlmfPbZZ0YvSxdr1qzB2LFjfc29ysrK8Oc//9noZemusrISgiBg8eLFRi8lLAYiBps6dSr+8Ic/4LPPPsMbb7yBhoYG3HrrrUYvK6b+/ve/w+v14sUXX8Snn36KVatW4Te/+Q0ee+wxo5emi97eXtx2221YsGCB0UuJqddeew2LFy/G448/jgMHDuC6667DrFmz0NTUZPTSYqq7uxtXXXUVVq9ebfRSdLVz504sXLgQe/bswdatW3Hu3DnMnDkT3d3dRi8t5oYPH45nn30W+/fvx/79+zFt2jTMnTsXn376qdFL082+ffvw0ksvYezYsUYvRZlIprJ582ZREASxt7fX6KXo6uc//7lYVFRk9DJ0tXbtWtHhcBi9jJiZOHGiOH/+/ICfXX755eKjjz5q0Ir0B0DcuHGj0cuIi6+++koEIO7cudPopcTF4MGDxd/+9rdGL0MXXV1d4ujRo8WtW7eK3/72t8WHHnrI6CWFxR0RE+no6MDvf/97fOtb30JqaqrRy9GVy+VCbm6u0csglXp7e/Hxxx9j5syZAT+fOXMmPvroI4NWRbHkcrkAIOn/v/R4PNiwYQO6u7tRVlZm9HJ0sXDhQsyePRszZswweimqMBAxgaVLlyIrKwt5eXloamrC5s2bjV6SrhoaGvDrX/8a8+fPN3oppFJbWxs8Hg+GDh0a8POhQ4eitbXVoFVRrIiiiCVLlmDKlCkoLS01ejm6OHToEAYOHAi73Y758+dj48aNKCkpMXpZMbdhwwbU1NSgsrLS6KWoxkBEB0899RQEQQj7tX//ft/tH3nkERw4cADvvvsuUlJS8P3vfx9iAjS81fo6AaC5uRkVFRW47bbb8MMf/tCglWsXyWtNRoIgBHwvimK/n1HiWbRoEQ4ePIhXX33V6KXo5rLLLkNtbS327NmDBQsW4O6770Z9fb3Ry4qpY8eO4aGHHkJVVRXS09ONXo5qbPGug7a2NrS1tYW9zciRI2XfKP/85z9RWFiIjz76yPTbhlpfZ3NzM6ZOnYpvfvObWLduHWy2xImDI/k7XbduHRYvXoyTJ0/qvDr99fb2IjMzE6+//jrmzZvn+/lDDz2E2tpa7Ny508DV6UcQBGzcuBE333yz0UvRzYMPPohNmzbh/fffR1FRkdHLiZsZM2aguLgYL774otFLiZlNmzZh3rx5SElJ8f3M4/FAEATYbDa43e6A35nFAKMXkIzy8/ORn58f0X2luNDtdsdySbrQ8jqPHz+OqVOnYsKECVi7dm1CBSFAdH+nySAtLQ0TJkzA1q1bAwKRrVu3Yu7cuQaujCIliiIefPBBbNy4Ee+9956lghCg7/Unwr+zWkyfPh2HDh0K+Nk999yDyy+/HEuXLjVlEAIwEDHU3r17sXfvXkyZMgWDBw/GP/7xDzz55JMoLi42/W6IFs3NzbjhhhswYsQI/PKXv8TXX3/t+53T6TRwZfpoampCR0cHmpqa4PF4UFtbCwAYNWoUBg4caOziorBkyRLcdddduOaaa1BWVoaXXnoJTU1NSZfrc+rUKRw5csT3fWNjI2pra5Gbm4sRI0YYuLLYWrhwIdavX4/NmzcjOzvbl+vjcDiQkZFh8Opi67HHHsOsWbNQWFiIrq4ubNiwAe+99x6qq6uNXlpMZWdn98vxkfIPTZ37Y1zBDh08eFCcOnWqmJubK9rtdnHkyJHi/PnzxX/+859GLy2m1q5dKwKQ/UpGd999t+xr3bFjh9FLi9rzzz8vXnLJJWJaWpo4fvz4pCz13LFjh+zf391332300mIq1P+Ta9euNXppMXfvvff63rcXXXSROH36dPHdd981ellxkQjlu8wRISIiIsMk1kE9ERERJRUGIkRERGQYBiJERERkGAYiREREZBgGIkRERGQYBiJERERkGAYiREREZBgGIkRERGQYBiJERERkGAYiREREZBgGIkRERGQYBiJERERkmP8PODHJSaQXmfUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(y,z)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regress y on x. The slope we get should be on average equal to $i = 0.9$. (See (29) in the Appendix section.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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.804</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.804</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>4.093e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Mon, 15 Jan 2024</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>20:41:33</td>     <th>  Log-Likelihood:    </th> <td> -5942.0</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td> 10000</td>      <th>  AIC:               </th> <td>1.189e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  9998</td>      <th>  BIC:               </th> <td>1.190e+04</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>const</th> <td>   -0.0002</td> <td>    0.004</td> <td>   -0.047</td> <td> 0.963</td> <td>   -0.009</td> <td>    0.008</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    0.8895</td> <td>    0.004</td> <td>  202.318</td> <td> 0.000</td> <td>    0.881</td> <td>    0.898</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.070</td> <th>  Durbin-Watson:     </th> <td>   2.015</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.586</td> <th>  Jarque-Bera (JB):  </th> <td>   1.089</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.006</td> <th>  Prob(JB):          </th> <td>   0.580</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.950</td> <th>  Cond. No.          </th> <td>    1.00</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/latex": [
       "\\begin{center}\n",
       "\\begin{tabular}{lclc}\n",
       "\\toprule\n",
       "\\textbf{Dep. Variable:}    &        y         & \\textbf{  R-squared:         } &     0.804   \\\\\n",
       "\\textbf{Model:}            &       OLS        & \\textbf{  Adj. R-squared:    } &     0.804   \\\\\n",
       "\\textbf{Method:}           &  Least Squares   & \\textbf{  F-statistic:       } & 4.093e+04   \\\\\n",
       "\\textbf{Date:}             & Mon, 15 Jan 2024 & \\textbf{  Prob (F-statistic):} &     0.00    \\\\\n",
       "\\textbf{Time:}             &     20:41:33     & \\textbf{  Log-Likelihood:    } &   -5942.0   \\\\\n",
       "\\textbf{No. Observations:} &       10000      & \\textbf{  AIC:               } & 1.189e+04   \\\\\n",
       "\\textbf{Df Residuals:}     &        9998      & \\textbf{  BIC:               } & 1.190e+04   \\\\\n",
       "\\textbf{Df Model:}         &           1      & \\textbf{                     } &             \\\\\n",
       "\\textbf{Covariance Type:}  &    nonrobust     & \\textbf{                     } &             \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lcccccc}\n",
       "               & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]}  \\\\\n",
       "\\midrule\n",
       "\\textbf{const} &      -0.0002  &        0.004     &    -0.047  &         0.963        &       -0.009    &        0.008     \\\\\n",
       "\\textbf{x1}    &       0.8895  &        0.004     &   202.318  &         0.000        &        0.881    &        0.898     \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lclc}\n",
       "\\textbf{Omnibus:}       &  1.070 & \\textbf{  Durbin-Watson:     } &    2.015  \\\\\n",
       "\\textbf{Prob(Omnibus):} &  0.586 & \\textbf{  Jarque-Bera (JB):  } &    1.089  \\\\\n",
       "\\textbf{Skew:}          &  0.006 & \\textbf{  Prob(JB):          } &    0.580  \\\\\n",
       "\\textbf{Kurtosis:}      &  2.950 & \\textbf{  Cond. No.          } &     1.00  \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "%\\caption{OLS Regression Results}\n",
       "\\end{center}\n",
       "\n",
       "Notes: \\newline\n",
       " [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.804\n",
       "Model:                            OLS   Adj. R-squared:                  0.804\n",
       "Method:                 Least Squares   F-statistic:                 4.093e+04\n",
       "Date:                Mon, 15 Jan 2024   Prob (F-statistic):               0.00\n",
       "Time:                        20:41:33   Log-Likelihood:                -5942.0\n",
       "No. Observations:               10000   AIC:                         1.189e+04\n",
       "Df Residuals:                    9998   BIC:                         1.190e+04\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         -0.0002      0.004     -0.047      0.963      -0.009       0.008\n",
       "x1             0.8895      0.004    202.318      0.000       0.881       0.898\n",
       "==============================================================================\n",
       "Omnibus:                        1.070   Durbin-Watson:                   2.015\n",
       "Prob(Omnibus):                  0.586   Jarque-Bera (JB):                1.089\n",
       "Skew:                           0.006   Prob(JB):                        0.580\n",
       "Kurtosis:                       2.950   Cond. No.                         1.00\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = sm.add_constant(x)\n",
    "model = sm.OLS(y, x).fit()\n",
    "predictions = model.predict(x) # make the predictions by the model\n",
    "# Print out the statistics\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, we regress y on x and z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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.849</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.849</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>2.809e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Mon, 15 Jan 2024</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>20:41:34</td>     <th>  Log-Likelihood:    </th> <td> -4632.3</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td> 10000</td>      <th>  AIC:               </th> <td>   9271.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  9997</td>      <th>  BIC:               </th> <td>   9292.</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>const</th> <td>    0.0020</td> <td>    0.004</td> <td>    0.520</td> <td> 0.603</td> <td>   -0.006</td> <td>    0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    0.4552</td> <td>    0.009</td> <td>   51.577</td> <td> 0.000</td> <td>    0.438</td> <td>    0.472</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>    <td>    0.4845</td> <td>    0.009</td> <td>   54.712</td> <td> 0.000</td> <td>    0.467</td> <td>    0.502</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 0.025</td> <th>  Durbin-Watson:     </th> <td>   1.995</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.988</td> <th>  Jarque-Bera (JB):  </th> <td>   0.031</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.004</td> <th>  Prob(JB):          </th> <td>   0.985</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.995</td> <th>  Cond. No.          </th> <td>    4.35</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/latex": [
       "\\begin{center}\n",
       "\\begin{tabular}{lclc}\n",
       "\\toprule\n",
       "\\textbf{Dep. Variable:}    &        y         & \\textbf{  R-squared:         } &     0.849   \\\\\n",
       "\\textbf{Model:}            &       OLS        & \\textbf{  Adj. R-squared:    } &     0.849   \\\\\n",
       "\\textbf{Method:}           &  Least Squares   & \\textbf{  F-statistic:       } & 2.809e+04   \\\\\n",
       "\\textbf{Date:}             & Mon, 15 Jan 2024 & \\textbf{  Prob (F-statistic):} &     0.00    \\\\\n",
       "\\textbf{Time:}             &     20:41:34     & \\textbf{  Log-Likelihood:    } &   -4632.3   \\\\\n",
       "\\textbf{No. Observations:} &       10000      & \\textbf{  AIC:               } &     9271.   \\\\\n",
       "\\textbf{Df Residuals:}     &        9997      & \\textbf{  BIC:               } &     9292.   \\\\\n",
       "\\textbf{Df Model:}         &           2      & \\textbf{                     } &             \\\\\n",
       "\\textbf{Covariance Type:}  &    nonrobust     & \\textbf{                     } &             \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lcccccc}\n",
       "               & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]}  \\\\\n",
       "\\midrule\n",
       "\\textbf{const} &       0.0020  &        0.004     &     0.520  &         0.603        &       -0.006    &        0.010     \\\\\n",
       "\\textbf{x1}    &       0.4552  &        0.009     &    51.577  &         0.000        &        0.438    &        0.472     \\\\\n",
       "\\textbf{x2}    &       0.4845  &        0.009     &    54.712  &         0.000        &        0.467    &        0.502     \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lclc}\n",
       "\\textbf{Omnibus:}       &  0.025 & \\textbf{  Durbin-Watson:     } &    1.995  \\\\\n",
       "\\textbf{Prob(Omnibus):} &  0.988 & \\textbf{  Jarque-Bera (JB):  } &    0.031  \\\\\n",
       "\\textbf{Skew:}          &  0.004 & \\textbf{  Prob(JB):          } &    0.985  \\\\\n",
       "\\textbf{Kurtosis:}      &  2.995 & \\textbf{  Cond. No.          } &     4.35  \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "%\\caption{OLS Regression Results}\n",
       "\\end{center}\n",
       "\n",
       "Notes: \\newline\n",
       " [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.849\n",
       "Model:                            OLS   Adj. R-squared:                  0.849\n",
       "Method:                 Least Squares   F-statistic:                 2.809e+04\n",
       "Date:                Mon, 15 Jan 2024   Prob (F-statistic):               0.00\n",
       "Time:                        20:41:34   Log-Likelihood:                -4632.3\n",
       "No. Observations:               10000   AIC:                             9271.\n",
       "Df Residuals:                    9997   BIC:                             9292.\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.0020      0.004      0.520      0.603      -0.006       0.010\n",
       "x1             0.4552      0.009     51.577      0.000       0.438       0.472\n",
       "x2             0.4845      0.009     54.712      0.000       0.467       0.502\n",
       "==============================================================================\n",
       "Omnibus:                        0.025   Durbin-Watson:                   1.995\n",
       "Prob(Omnibus):                  0.988   Jarque-Bera (JB):                0.031\n",
       "Skew:                           0.004   Prob(JB):                        0.985\n",
       "Kurtosis:                       2.995   Cond. No.                         4.35\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = sample[:,[0,2]]\n",
    "X = sm.add_constant(X)\n",
    "model = sm.OLS(y, X).fit()\n",
    "# Print out the statistics\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficients for x and z should be on average equal to $\\frac{i-jk}{1-j^2}$ and $\\frac{-ij + k}{1-j^2}$. (See (21) in the Appendix.) In our case, $\\frac{i-jk}{1-j^2} = \\frac{-ij + k}{1-j^2} = 0.473$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.47368421052631576"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.9 - 0.9*0.9)/(1-0.9**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, in this particular case, including Z in the regression indeed makes the coefficient on X smaller. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if this finding can be generalized by looping over all posible values of $i$, $j$, and $k$ from -0.9 to 0.9 with 0.1 increments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = [i * 0.1 for i in range(-9,10)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "OUTPUT = []\n",
    "for i in values:\n",
    "    for j in values:\n",
    "        for k in values:\n",
    "            mean = np.array([0, 0, 0])\n",
    "            cov = np.array([[1, i, j], [i, 1, k], [j, k, 1]])\n",
    "            eigenvalues = np.linalg.eig(cov)[0]\n",
    "            if (eigenvalues[0] >= 0) and (eigenvalues[1] >=0) and (eigenvalues[2] >= 0):\n",
    "                ps = True\n",
    "                sample = np.random.multivariate_normal(mean, cov, 10000)\n",
    "                x = sample[:,0]\n",
    "                y = sample[:,1]\n",
    "                z = sample[:,2]\n",
    "\n",
    "                x = sm.add_constant(x)\n",
    "                model = sm.OLS(y, x).fit()\n",
    "                b1 = model.params[1]\n",
    "\n",
    "                X = sample[:,[0,2]]\n",
    "                X = sm.add_constant(X)\n",
    "                model = sm.OLS(y, X).fit()\n",
    "                b2 = model.params[1]\n",
    "                beta = (i - j*k)/(1 - j**2)\n",
    "                gamma = (-i*j + k)/(1 - j**2)\n",
    "\n",
    "                if b1 > b2:\n",
    "                    est = \"overestimated\"\n",
    "                if b1 < b2:\n",
    "                    est = \"underestimated\"\n",
    "            else:\n",
    "                ps = False\n",
    "                b1 = None\n",
    "                b2 = None\n",
    "                est = None\n",
    "                beta = None\n",
    "                gamma = None\n",
    "\n",
    "            row = [i, j, k, b1, b2, est, beta, ps, gamma]\n",
    "            OUTPUT.append(row)         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And save the results to a dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = pd.DataFrame(OUTPUT, columns = [\"i\",\"j\",\"k\",\"b_xy\", \"b_xyz\", \"estimate\", \"beta\", \"positive semidefinite\", \"gamma\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show only the rows with a positive semidefinite covariance matrix. \"b_xy\" is the estimated regression slope for x without controlling for z. \"b_xyz\" is the estimated regression slope for x after controlling for z. $\\beta = \\frac{i-jk}{1-j^2}$  and $\\gamma = \\frac{-ij + k}{1-j^2}$ are the theoretical slopes for X and Z. You can see from this table that \"b_xy\" is a very close estimate of $i$, and \"b_xyz\" is a very close estimate of $\\beta$. \n",
    "\n",
    "Without controling for z, bias = $\\gamma j$. When $\\gamma j$ is postive, \"b_xy\" overestimates $\\beta$. When $\\gamma j$ is negative, \"b_xy\" underestimates $\\beta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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>i</th>\n",
       "      <th>j</th>\n",
       "      <th>k</th>\n",
       "      <th>b_xy</th>\n",
       "      <th>b_xyz</th>\n",
       "      <th>estimate</th>\n",
       "      <th>beta</th>\n",
       "      <th>positive semidefinite</th>\n",
       "      <th>gamma</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.9</td>\n",
       "      <td>0.7</td>\n",
       "      <td>-0.902897</td>\n",
       "      <td>-1.410553</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>-1.421053</td>\n",
       "      <td>True</td>\n",
       "      <td>-0.578947</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.9</td>\n",
       "      <td>0.8</td>\n",
       "      <td>-0.895188</td>\n",
       "      <td>-0.926884</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>-0.947368</td>\n",
       "      <td>True</td>\n",
       "      <td>-0.052632</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>-0.903159</td>\n",
       "      <td>-0.478150</td>\n",
       "      <td>underestimated</td>\n",
       "      <td>-0.473684</td>\n",
       "      <td>True</td>\n",
       "      <td>0.473684</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.8</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.899779</td>\n",
       "      <td>-1.389657</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>-1.388889</td>\n",
       "      <td>True</td>\n",
       "      <td>-0.611111</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.8</td>\n",
       "      <td>0.6</td>\n",
       "      <td>-0.900013</td>\n",
       "      <td>-1.170034</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>-1.166667</td>\n",
       "      <td>True</td>\n",
       "      <td>-0.333333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6838</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.8</td>\n",
       "      <td>0.8</td>\n",
       "      <td>0.909854</td>\n",
       "      <td>0.725262</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>0.722222</td>\n",
       "      <td>True</td>\n",
       "      <td>0.222222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6839</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.8</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.904283</td>\n",
       "      <td>0.506119</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>True</td>\n",
       "      <td>0.500000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6856</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.7</td>\n",
       "      <td>0.900043</td>\n",
       "      <td>1.420639</td>\n",
       "      <td>underestimated</td>\n",
       "      <td>1.421053</td>\n",
       "      <td>True</td>\n",
       "      <td>-0.578947</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6857</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.8</td>\n",
       "      <td>0.896327</td>\n",
       "      <td>0.952184</td>\n",
       "      <td>underestimated</td>\n",
       "      <td>0.947368</td>\n",
       "      <td>True</td>\n",
       "      <td>-0.052632</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6858</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.899231</td>\n",
       "      <td>0.462386</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>0.473684</td>\n",
       "      <td>True</td>\n",
       "      <td>0.473684</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>4775 rows × 9 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        i    j    k      b_xy     b_xyz        estimate      beta  \\\n",
       "16   -0.9 -0.9  0.7 -0.902897 -1.410553   overestimated -1.421053   \n",
       "17   -0.9 -0.9  0.8 -0.895188 -0.926884   overestimated -0.947368   \n",
       "18   -0.9 -0.9  0.9 -0.903159 -0.478150  underestimated -0.473684   \n",
       "33   -0.9 -0.8  0.5 -0.899779 -1.389657   overestimated -1.388889   \n",
       "34   -0.9 -0.8  0.6 -0.900013 -1.170034   overestimated -1.166667   \n",
       "...   ...  ...  ...       ...       ...             ...       ...   \n",
       "6838  0.9  0.8  0.8  0.909854  0.725262   overestimated  0.722222   \n",
       "6839  0.9  0.8  0.9  0.904283  0.506119   overestimated  0.500000   \n",
       "6856  0.9  0.9  0.7  0.900043  1.420639  underestimated  1.421053   \n",
       "6857  0.9  0.9  0.8  0.896327  0.952184  underestimated  0.947368   \n",
       "6858  0.9  0.9  0.9  0.899231  0.462386   overestimated  0.473684   \n",
       "\n",
       "      positive semidefinite     gamma  \n",
       "16                     True -0.578947  \n",
       "17                     True -0.052632  \n",
       "18                     True  0.473684  \n",
       "33                     True -0.611111  \n",
       "34                     True -0.333333  \n",
       "...                     ...       ...  \n",
       "6838                   True  0.222222  \n",
       "6839                   True  0.500000  \n",
       "6856                   True -0.578947  \n",
       "6857                   True -0.052632  \n",
       "6858                   True  0.473684  \n",
       "\n",
       "[4775 rows x 9 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results[(results['positive semidefinite'] == True)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some interesting cases:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When X and Y, Y and Z are correlateed at 0.9, varying the correlation betweeen X and Z from 0.9 to 0.7 would vastly change the result of control.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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>i</th>\n",
       "      <th>j</th>\n",
       "      <th>k</th>\n",
       "      <th>b_xy</th>\n",
       "      <th>b_xyz</th>\n",
       "      <th>estimate</th>\n",
       "      <th>beta</th>\n",
       "      <th>positive semidefinite</th>\n",
       "      <th>gamma</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>6856</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.7</td>\n",
       "      <td>0.900043</td>\n",
       "      <td>1.420639</td>\n",
       "      <td>underestimated</td>\n",
       "      <td>1.421053</td>\n",
       "      <td>True</td>\n",
       "      <td>-0.578947</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6857</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.8</td>\n",
       "      <td>0.896327</td>\n",
       "      <td>0.952184</td>\n",
       "      <td>underestimated</td>\n",
       "      <td>0.947368</td>\n",
       "      <td>True</td>\n",
       "      <td>-0.052632</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6858</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.899231</td>\n",
       "      <td>0.462386</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>0.473684</td>\n",
       "      <td>True</td>\n",
       "      <td>0.473684</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        i    j    k      b_xy     b_xyz        estimate      beta  \\\n",
       "6856  0.9  0.9  0.7  0.900043  1.420639  underestimated  1.421053   \n",
       "6857  0.9  0.9  0.8  0.896327  0.952184  underestimated  0.947368   \n",
       "6858  0.9  0.9  0.9  0.899231  0.462386   overestimated  0.473684   \n",
       "\n",
       "      positive semidefinite     gamma  \n",
       "6856                   True -0.578947  \n",
       "6857                   True -0.052632  \n",
       "6858                   True  0.473684  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results[(results['positive semidefinite'] == True)].tail(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some cases, controling for Z wouldn't matter at all even though Z is correlated with both X and Y. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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>i</th>\n",
       "      <th>j</th>\n",
       "      <th>k</th>\n",
       "      <th>b_xy</th>\n",
       "      <th>b_xyz</th>\n",
       "      <th>estimate</th>\n",
       "      <th>beta</th>\n",
       "      <th>positive semidefinite</th>\n",
       "      <th>gamma</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>450</th>\n",
       "      <td>-0.8</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.4</td>\n",
       "      <td>-0.803613</td>\n",
       "      <td>-0.803792</td>\n",
       "      <td>overestimated</td>\n",
       "      <td>-0.8</td>\n",
       "      <td>True</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       i    j    k      b_xy     b_xyz       estimate  beta  \\\n",
       "450 -0.8 -0.5  0.4 -0.803613 -0.803792  overestimated  -0.8   \n",
       "\n",
       "     positive semidefinite  gamma  \n",
       "450                   True    0.0  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results[(results['i'] == -0.8) & (results['j'] == -0.5) & (results['k'] == 0.4)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Z $\\nsim$ X is no guarentee for Z not being a confounder for the relationship between X and Y. This is probably the most counter-intuitive part of the results. "
   ]
  },
  {
   "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>i</th>\n",
       "      <th>j</th>\n",
       "      <th>k</th>\n",
       "      <th>b_xy</th>\n",
       "      <th>b_xyz</th>\n",
       "      <th>estimate</th>\n",
       "      <th>beta</th>\n",
       "      <th>positive semidefinite</th>\n",
       "      <th>gamma</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.9</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.8</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.7</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.6</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>85</th>\n",
       "      <td>-0.9</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6773</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6792</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.6</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6811</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.7</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6830</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.8</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6849</th>\n",
       "      <td>0.9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>361 rows × 9 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        i    j    k  b_xy  b_xyz estimate  beta  positive semidefinite  gamma\n",
       "9    -0.9 -0.9  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "28   -0.9 -0.8  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "47   -0.9 -0.7  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "66   -0.9 -0.6  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "85   -0.9 -0.5  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "...   ...  ...  ...   ...    ...      ...   ...                    ...    ...\n",
       "6773  0.9  0.5  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "6792  0.9  0.6  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "6811  0.9  0.7  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "6830  0.9  0.8  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "6849  0.9  0.9  0.0   NaN    NaN     None   NaN                  False    NaN\n",
       "\n",
       "[361 rows x 9 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results[(results['k'] == 0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lessons\n",
    "\n",
    "Bivariate associations can tell neither whether a third variable is a confounder nor the direction of bias it causes. What determines these issues are conditional associations rather than unconditional associations. See the appendix for a mathematical proof. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<img src=\"Confounder_Page_1.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<img src=\"Confounder_Page_2.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}