week09/Week09.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate regression\n",
    "\n",
    "Up until now we have explored simple linear regression (SLR).\n",
    "SLR supposed a the conditional distribution of a random variable Y has a Normal distribution where the mean is a function of a random variable X and the variance is constant.\n",
    "\n",
    "But often there may be many variables that change the probability distribution of Y.\n",
    "\n",
    "Multivariate Linear Regression (MLR) is a statistical model that relates more than one random variable X (called covariates) to Y.\n",
    "\n",
    "We will first rewrite SLR using a centered $X$ random variable $X*$ to simplify maximum likelihood estimates of $\\beta_{0}$, $\\beta_{1}$. Next, we will rewrite these MLEs using matrix algebra. Finally, we generalize these MLE estimating equations to more than one variable $X$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Centering on X\n",
    "\n",
    "Our typical SLR setup assumes a dataset of $x,y$ pairs generated from pairs of random variables and a specific probability distribution for $Y$. Suppose we collect a sample of $N$ data points $\\mathcal{D} = [(x_{1},y_{1}),(x_{2},y_{2}),\\cdots,(x_{n},y_{n})]$ which are generated from $N$ pairs of random variables $(X,Y)$.\n",
    "\n",
    "SLR supposes the following conditional probability distriution for $Y$\n",
    "\n",
    "\\begin{align}\n",
    "    Y|x,\\beta,\\sigma^{2} \\sim \\mathcal{N}\\left(\\beta_{0}+\\beta_{1}x,\\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "or\n",
    "\\begin{align}\n",
    "    Y_{i} =  \\beta_{0}+\\beta_{1}x_{i} + \\epsilon_{i}\\\\\n",
    "    \\epsilon_{i} \\sim \\mathcal{N}\\left(0,\\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "\n",
    "Lets change the above SLR model by creating a new variable $x*$by subtracting from each $x_{i}$ the sample average $\\bar{x}$.\n",
    "This is called centering on x and $x^{*}_{i} = x_{i} - \\bar{x}$ is called a centered covariate.\n",
    "\n",
    "The SLR model relating $x^{*}$ and $Y$ is then \n",
    "\n",
    "\\begin{align}\n",
    "    Y_{i} =  \\beta_{0}+\\beta_{1}x^{*}_{i} + \\epsilon_{i}\\\\\n",
    "    \\epsilon_{i} \\sim \\mathcal{N}\\left(0,\\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "\n",
    "We can visualize the effect of centering below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x1c902af70>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAEjCAYAAACmQt93AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4I0lEQVR4nO3de3TT9f0/8OfnkqRN7y1tQKBlYJHbEMVBdVM2VNQ5FO+DKYrzMnD7ThSP6Jy7fDcFRZnfTdm+MuZh07PxVUBlHHacMqeTy76/LxflWoWCILS0TW9Jm8vn8/79kSYkaXpJmzTJ5/N8nMM5NE2Td5s++8rnfZWampoEiIiITEBOdQOIiIgGC4seERGZBoseERGZBoseERGZBoseERGZBoseERGZBoseERGZBoseERGZhprqBlBiFf7hZMTHTQuGp6glqVddXY3KyspUN4PiZNbXjdk9K5m/A7zSIyIi02DRIyIi02DRIyIi02DRIyIi02DRIyIi0+DsTaIUqXP7sbbajVq3BoddwfxKO8rsjCRRdxKRGSaMKAXq3H4s3tYMjy6gSBKOtmnYVe/DyosLWPiIYkhUZti9SZQCL+1rw0GnF4edPhxp9kHTBTy6wNpqd6qbRpSWEpUZFj2iQVbn9mNDTQeafUC7Bji9Avsa/dB0gbp2LdXNI0o7icwMix7RIFtb7YYqAxACACAB0AEcb/WjLFtJZdOI0lIiM8PBA6JBEhyEX3/EDU3TAQkQQkCSJEAIaELG/Ep7qptJlDaSkRkWPaIE6m52WfggvFcHWv0AhEC+TYYmJFhkCbPLbZzEQqYTKzPB25ORGSaMKEF6ml22ttodun14joIWr4AGQJUkjClQYZMlLJqYl+pvgWhQdZeZHziAN5OUGY7pESVIeGEDAEWSQrPLat1a6PYsVcboPBkCQINHg08T+NGUXF7lkel0l5mNtWrSMsOiR5Qg4SENUiQJde2BbhutcxC+w6/jSKsOCUCJTYFFkfDL3W2oc/tT0Gqi1OkuM40+KWmZYdEjSpDwkAZpQqAsOzBOYZMlaELgpEuDJgLvbkfkKhFXhERm0l1mii0iaZlh0SNKkPCQAoHw2mQpNJll5cUFqCqzwqYAxTYZE4tVWJVABINXhERm0l1m5jj8ScsMBxGIEiQY0rXVbtS1a6ErvOC4Q5ldxZLz8wEA2+u8Ed06wStCIjPpLjPNJxtDn090Zlj0iBIoPKTdmV9px656X2gAP/yKkMhsYmWmOeo+icwMix7RIOvtipCIIiUyM0wZURJ1t1i9L1eERGZU5/bj98dVaPXOpGSGRY+oBwM5v4vHB5EZJSIzzhYFeUJPSmY4e5OoG8EAbq/z4rhLx/Y6LxZva+7z2qCeFqsTGVGiMiN3zldJRmZY9Ii6MdCi1dNidSIjyoTMsOgRdWOgAexpsTqREWVCZlj0iLox0ADGWngrdIE2r45HtjmxYk8Ltx4jQ0lUZnRx9msTnRkWPaJu9LTDSl+E7yhRkStjUqECWQI+afL3a7yDKN0lKjOT87WkZYZFj6gb0UWrqswa9yyy4DTrZ6qKkGtVICSJE1vIsBKVme+W+5OWGc6bJupBItfTcWILmUG6Z4ZXekSDhBNbiOKTjMzwSo8oQXpblMs9N4kipSIzLHpECdCX3Ve45ybRWT1lJigZmWHaiBKgp0W54eMb3HOTKKCnzFwfdiGX6MxwTI8oAThJhSg+qcoMix5RAnCSClF8UpUZFj2iBBjoolwis0lVZjimR5QAnKRCFJ+eMhN9cnoiMZFECcJJKkTxSUVm2L1JRESmwaJHRESmwaJHRESmwTE9ojC9bYtERJEyLTPp2zKiQdaXrcSI6KxMzAy7N4k69bQtEhF1lYmZYdEj6sStxIjik4mZSc/rT6JBEj4ecbjZDwGBbPXsNkjcSowoUqZnJmOu9J5//nl84xvfwMiRIzFmzBjcdttt2L9/f8R9hBB4+umnMW7cOAwdOhTXXnstDhw4kKIWUzqoc/uxYk8LHtnmxIo9Lahz+yM+t3hbM7bXeXHcpUMXAoeaNLT7A+9SM30rMWaG+sPomcmYovfhhx/iu9/9Lv72t7/hrbfegqqqmDNnDpxOZ+g+L7zwAl588UUsX74c7733HkpLS3HDDTegtbU1hS2nVGnwICKg2+u8WLytORTi6PEIu0XB2AIFEiRU5MqoKrOm9YB8b5gZild0UTNiZtK3ZVHWr18f8fHvfvc7lJeXY/v27bjmmmsghMCqVavw4IMP4vrrrwcArFq1CpWVlXj99dexYMGCVDSbUmhjrdrjGXexxiPsFgUVuTKeqSpKRZMTipmhePV2LqQRMpMxV3rR2traoOs6CgsLAQDHjh1DbW0tZs6cGbpPdnY2LrnkEuzYsSNFraRUavBJPQ6ym+04IGaGetPbxBQjZCZjrvSiLV26FF/+8pcxbdo0AEBtbS0AoLS0NOJ+paWlOHXqVLePU11dnbxGpkRkX7rxvr++K7GoONnighyWYV0AsqShuroeX1OAD91WuHUJshT4nFUW+JriRHV1beoaHqaysjJhj5VJmTHn723qs6u4VbS2KmmRmf5+/71lJiOL3uOPP47t27djy5YtUJSBvcNI5B+VtPDhyYgPDff9xWGOpxqfIz/UXRMcZH+oc8yhEsB/j/ab4jigTMpMdXW1OX9v0yC7Dw2PXGyeqswk83cg49L92GOPYf369Xj77bcxatSo0O0OhwMAcObMGYwcOTJ0+5kzZ1BWVjbYzaQ0UGJDr2fcdXe0SaZtrdQTZob6qi/nQmZ6ZtKvRT149NFHsWHDBrz99tsYO3ZsxOcqKirgcDiwdetWXHjhhQCAjo4ObNu2DT//+c9T0VxKA/05rysTt1bqDjND8TJ6ZjJmIsuSJUvw2muv4eWXX0ZhYSFqa2tRW1uLtrY2AIAkSVi4cCFeeOEFvPXWW9i/fz8WLVqEnJwc3HzzzSluPWWSTNxaKRZmhgZLJmUmvUpwD1avXg0AoanVQY8++igee+wxAMAPf/hDtLe345FHHkFTUxOmTp2K9evXIy8vb9DbS5krE7dWioWZocGSSZnJmKLX1NTU630kScJjjz0WCjRRfzjsCo62RYY406ZlA8wMDZ5MykzGFD2iZIg1+D6/0o5d9b4uM9jSeWslosGS6Zlh0SNTqXP78dK+Vnx42guvDrT5dAzPUWC3KBGD773NYCMyC6NlJv1aRJQkdW4/Fn7gxJ5GHyAAt1/AqwPNXoFJxRKyVDliy6V4Z7ARGY0RM5MxszeJBmpttRs1bRogAmNZOgBJCgT5pCsw4J6ug+9EqWDEzLDokWnUujV49UB4AUDpHHMXALx64P/pOvhOlApGzAyLHpmGw67AKgfOkAOALEWCBEACYJUz4ywwosFkxMyw6JFpzK+0Y1SuAkiBEEuShDwLcI5dQVWZJSPOAiMaTEbMTOa0lKifwqdYf7nYgtF5MnY1+AFJwtcc2Vg0MTejQkuUbEbOTGa2mgh92+A2ek/AYHfMa5eXZGxoifqLmWH3JmWoYDC313lx3KVje50Xi7c1o87tj7hfJu0JSJRMzEwAix5lpL4GM5P2BCRKJmYmIPOvVcmUwoPZ4ddx0qXBowNNHVpEl00m7QlIlEzRxSyYm1OuwJVeMDdGzwyv9CgjOewKNCHQ4ddxwOlHk0dHu1/A6dEjumzmV9phkwPjEkBmTrEmSoRgZgCEcuP06PBqiOjqNHpmeKVHKdefE5eDG9xWO71o9QkIAAoExuaroS6b6+19OwmaKNMMJDMtXg37nH64fAKqBJQWSBFdnUvOzzd0ZozxXVDG6u3E5e7CXWZX8aMpufjmFg80HfAH3pTi32f8OL9EoK5dATrfmPbnJGiidDXQzNz1vhPtPgGfDvhwNjMl2ZbQuJ2RM8PuTUqpngbXe5tttvmEB3kq4BOBbZEEAv//f/Ua0Nk1Q2Q0A81MiU2CNyoz/1evocXjN8y4XU9Y9Cilepop1ttss1q3htbI2dYIPtLeBl+ym06UEgPNzKFmDeFfLSFQ/A43a4YZt+sJuzcppXqaKdZTuOvcfhxu9qPNKyICLABkyYDTI/D74yq0emefxzyIMsFAM9PiEdClwObReucVnyIBBSriHifMRLzSo5TqaaZY+GyzIE0I2BUJCz9woqbFBw2AjkBwZQmwyIF/9R4de1uUHhfhEmWihGRGAJo4mxm7AjT70evCdSNg0aOUCs6urCqzoiJXjtjAtrtwu3wa9jT64PIDeZ1vRAUCw3h2BdAAnJsvQ+58w2u0HSXI3JKZGaPuwhLOeNeulHG6mykWa7nBN0fYcPM7DWj1Cqhy4KiTIhlo9QlYZOD8IVYMsUpo0yW42ryhxzLSjhJEycpMOKNmhkWP0lp4uIMz09r8AjoCh1j6dIE8i4RCm4xReQpevqwI9/7Tic9afJA1CV/K0pGlyobaUYKoJwPJTJYiYXiOYujMsOhRyvV1oW1wZlqxTcZJlx444wtAu18gRwUuLFGxeFszBAS8GuD1STjg9OO8QgX5VsUUM9PIHJKVmXa/jhavMHRmWPQopXpbaBsuODOtIk9Fq88Htz+wE4sqAVOGWGBXJXh0HdmqgonFEo40+qErArIkZdxBl0TdSWZmTrRp6NB0Q2eGE1koZercftz7Tyf2NHhR0+JHh1/vcQA9ODMtS5UxqdiC4TkKhmTJuKDUihe/VgSXH6GBeKsiY2S2wIRiG84rVA0ZXjKfZGdmdIHF8Jlh0aOUCL5bPdLsg0cDmjyBDXCDIY41gB6cmeb2aTjp0tChBaZiP/2V/NDu8LGmaxtxXILMh5lJDBY9SongWIOqyHD5dLT5Bdp8Oo63+bsNXXDvwNp2HR5NR7YiYWiOjF/ubou5O7wuYKjd4cncBiMzRjtRIRYWPUqJWrcGTRdo9erwdS6U9QnglEuHJES3odt8woPyPBUTim0YXWBBtqqEunai1y9NztcMOy5B5jMYmQlf82dUxv3OKK057AreO9kBSBLyLECHJqDpQLZFwsSi7scTejvVOXy6dnV1vaHDS+YyGJkxA/5FoKSLNb16fqUdfzzsBkRgPMKuALIqYXyRCncP62HD9x0MnvzcoQl0+ANHqrDIkVGs2NPCzCSBub5bGnSxpldvO+3Bl4stKLAALT4gWxawqzKG5yiwKIEwRwc+GMzwgzAPNWnQROBxBQQWb2s2fNcMmcf2Oi8zkwQc06Okij7qxKcJ7Gn04e1j7SjJVqAisEtEMLxCF9jn9He78W1wDEKWJNiUwKLbicVqxDgFkREwM8nBokdJFT2ecNKlAZ0D8ACQa5Xh9QtUN/kwqVDBl4stEJLU48a3ZXYVYwvODsxbFTl0XyPuFUjmxswkFoseJVX0OiCPDqAzoPsa/Wj1CdhUGbIsobpFR32H3uOge3ePCxh/fRGZEzOTWCx6lFTR64CsMiAJASCwAa4EQAiBLEVCi1fDh6c92NfoxWfNPnT4dQCxg2nG9UVkLsxMcphn9JJSIvqok/EFNuxz+nHQ6YMECUIE9sYcYpNwqEmDBB1tfsCvA6fcGiYWKSjNtnQJZqwjVIx60jOZU1WZlZlJAnN9t5QS0euA6tx+zH+vAfubAgPtxTYJp9t1+HQdXh2wKRIAAb8OfNai41cX58YMptnWF5G5MDPJwe5NGlR1bj9e2teKE24dfh1QZQltPoHT7To6tEB4FUmCXZWRb5WRrUrYfMKT6mYTpQwzk1i80qNBE1yzd7DJhw6/QJYqwaMJFFoDB1b6dUQMyAfGLWTTzS4jCmJmEo9XejRogmv2vDogdc5Gy1Zl2FUJk0usUGQJeucge3DcYphdNt3sMqIgZibxWPRoUNS5/dhyvAOHnT64fDo0PTDLTALg0QUsioRbvpSFEpsMmwIU2mSMypNQ69ZR0+LHij0tocW2RGbAzCQHuzcp4aL32vzmCBt+ubsNTq8eOLlZCLj8QK6qQ5ZlWKTAEUCPTMnHI1MC725rWnzYecaHoTkyzngETtd5uz0dmijTxSpOi7c1MzNJwCs9SqjgGET4lkh3/cOJVp+GETkK5M4umhwVEAByVAmzK7JDwQzOLhuVb0F5XmCrJCD2LhNERhDMTDRmJjlY/imhovfaVCQJbk3glEvH6AILxhepOOnS4NWBsiwJr11eEvNdaG/HoRAZRTAz0ZiZ5GDRo4SKFbwsRUKHFhiPyFJljCkIzDyrKrN22+0SPA5F0wVOtGnwaAIWGRhfkJ3074FoMMXKDABmJknYvUkJlaMCnzb7sN/pC22LNMwuI1uV49r+aH6lHZIQ+KTBB6dXwO0XaPYJfNzo4+A8GUowM9GYmeRg0aOEqXP7sc/pR4tHR7tfoMmjY5/TB6ss4ZUZRagqs6IiV0ZVmbXXwfUyu4qJRSrybYHp2YU2GROLLJBkiWMUZBjhmYnGzCRHRhW9f/3rX/j2t7+N8ePHo7CwEK+++mrE54UQePrppzFu3DgMHToU1157LQ4cOJCi1prP2mo3hCRhUokFRVYJdouEHFWC26/jD4cDoVsyOQ9Lzs/v02wylx84t8CC8UUWjCmwIEuVOUYRJ2YmvYVnJhozkxwZVfRcLhcmTJiAZcuWITu7az/1Cy+8gBdffBHLly/He++9h9LSUtxwww1obW1NQWvNJzg2YVVkjC6wYEy+inY/8LlLj3m4ZW94FMrAMTPpLTwz0ZiZ5Mioojdr1iw8+eSTuP766yHLkU0XQmDVqlV48MEHcf3112PChAlYtWoV2tra8Prrr6eoxeYSHbgTbRq0ziNQgPinUPMolIFjZtJbrCIVjplJvIwqej05duwYamtrMXPmzNBt2dnZuOSSS7Bjx44Utsw8ogPn8uno0ARcfj00qSWerpbgUSjxjGtQ3zEzqRedmXDMTHIY5idRW1sLACgtLY24vbS0FKdOner266qrq5ParsEX+Y5usL+/HziAjbUqvuiQ4PYpkDUJHT6Bdq+ORrcPo7MFviRpqK6u7/NjXm9H6NtqPgl0XcbbvUx/fSsrK5P22OmcmUx/3eIRzMzfYY24vd7tT0lm0kV/fwd6y4xhil5/JfOPSkp8eDLiw8H+/ioBVE0CVuxpQZPSgYNOLXTasy4EWiQZD108ZFDeeVZXVxvv9U0Dyf6Zmu11C2bmt3+IzK5FVQY9M+kimb8DcXVvvvvuuxA99D+nksPhAACcOXMm4vYzZ86grKwsFU0yrcBGue042qojxwLkKgLZKlBskzGth8W1RsTMUF/Emqhi1swkW1xF7+abb8aECRPw4x//GB9//HGy2tQvFRUVcDgc2Lp1a+i2jo4ObNu2DdOnT09hy8wluI+gs3OtXqtXoNUPSEKgXRP43KWZaqEsM0O96W7vTbNmJtniKnqvvvoqpk+fjtWrV2PGjBm45JJL8Otf/7rH/v9Eamtrw969e7F3717ouo4TJ05g7969+PzzzyFJEhYuXIgXXngBb731Fvbv349FixYhJycHN99886C0z+jq3IHjSh7Z5uz22JLgPoLleSpkADqAFq9AXYeAVwMERFxTsDMdM2Nu8WSmy9eaNDPJJjU1NcXd99La2oqNGzfif/7nf/Cvf/0LkiTh0ksvxbe//W3Mnj0bdntypsd+8MEHmD17dpfb586di1WrVkEIgWXLluGVV15BU1MTpk6dihUrVmDChAlJaU86KowaF2haMDwhjxt8NxrcTDo4FTo4Myx4nND6I274dGBEbmBd0J4GH9x+AZsMXFhqhVU5u4fgkvPzE9K27qTT2BAz03fp9LoNRLyZqW6JnKFZYBn8zKSLZP4O9Kvohfviiy/w+uuvY926ddi/fz/sdjuuvfZazJs3DzNmzEhUO6mPklX0VuxpwfY6b8TGuMEgzq+0h8Jd0+KH0xOYZj2xWMVnzX64/QKFNhljCs7uOlGRK+OZqqKEtK076frHk5npWbq+bvGKNzMNnsg/xefmK4OemXSRNhNZYtE0DT6fD16vF0IIZGVl4f3338ecOXNw6aWXYv/+/YloJ6VYcOcIr6bjSLMPBxq9ONbiR02LP+I4oeE5Suhd7Yk2DRYZgAQMzzm7I4TZd4hgZswh3sxEY2aSo19Tgpqbm7Fx40b85S9/wY4dO6CqKmbNmoWf/OQnuOqqqyDLMjZv3ozHH38cDzzwQMRAOWWG6NPPc1QJ7X4tYgmCy69jZ13gD3cwtFmqHDr/y6oAVw7PxseNPkhy4PNm3SGCmTG+gWam7pQ34vEsirkzkyxxFb1NmzZh3bp1eOedd9DR0YELL7wQy5Ytw80334yiosjL7m9961tobGzEww8/nNAGU/JFj0UcbdMgdIGTbTo0ISBLUii0DruME24dFkWKCPGofCk0BhH8Y1DXrqEsW8H8SrtppmAzM+aQiMxEqyqzmjIzyRbXT/GOO+7AOeecg+9973uYO3cuxo4d2+P9J06ciFtuuWVADaTBF336uaYLHG8NjDtAkmCVBbJVGcNzFGSpMvIUYG+jD24tsM/mMLuMfKsSemdaZldNMwAfjZkxh0RkJppZM5NscRW9DRs2YMaMGZBi9D/HMnXqVEydOrVfDaPUCC4sr+0IzLgsy5LwaYsOHYBfSFAg0OQFAIEvXBqGZOk42iEwNEfGKZeODk1HbTuwfBr3+wOYGTNIVGZocMT1V+nrX/96kppB6SByYTnQLgROuQWsMgAR+NfkC/y33a/DIgM1rcDEIgnZqgWjCwID7ZoQWHfEjdwTntD4hlm7Z5gZY0tkZg43Ry5ZWLGnxbS5SSb+NCkkfGH5vkY/dEmCXxfQNECXAlN9FQnwC4QG5iUA1S0CBTY9NC7h0wQ21nhwbqEaGt/YVe/Dj6bkYjMLIRlIojITa/bm9jovdtR6MLFIhcsPZiZBDHO0EA1ccIo1AORYJPi0zjO5AOSoACRAF2eDK0uARQH8OnDSdfZd6gmXBkU6O76hSBJafRru+ocT2+u8/TockygdJSozsY4W0nSB3fU+vH3cw8wkEIsehTjsCtr9GvY1+tHqE7AqEuxqoGtGkiQoYW9GVTnwztWmSFBloEM7e2ilXwfK8yLfjZ5y6XBrkYUwnsMxidJRojJjk7te6Z1o0yAkCT498DEzkxgsehQyv9KO0+7AFGsJgZO1LbIMR3bg/4VWGYoMqBKAznevqiRhcrGKMflq6NDKG0ZlQYkKcYdfD52gHhTP4ZhE6ShRmVl5cUGXx/ZoAhCd44OdmJmBY+cwhZTZVXxliAU7zvjg1QGrLIV2hTjWquHcQhWjfDp21fugA3BkSSjPV5BnUUKhXVvtxpl2Dcdb/Rhql5GtKtCEQLZFhiM78j0Wd5mgTJeozDy7p7XLY3M3o+Rg0aMIXyqwoNYjEL1f4JxRNuRaFdS1a/hKqRWAgFtDaOEsgIjFuY5sGaddOqaVqRiVb8HPL8zDL3e3ddl8l7tMUKZLVGaiza7gbkbJwKJnItHbJMWaCTa/0o5d9b4uxWnRxLweZ42t2NMSEV67RUF5voxR+WcXpq+8WDHtziyUmQYzM9F++pVCU+9mlCz86ZlE+DZJmi7w3skO/PGwGzeMskWEs8yuYuXFBXhpXys+PO0FIOFrQ629Pn74LLag6PEHM+/MQpknFZmJxswkHieymERwPZGmC+xr9KPZB7h8Ot4+7ok5Dbq6RUdhloKSbAWfNPl7nSrtsCtdpl1z/IEyWSoyQ8nHomcSwXeVJ9rO7vgudU6Hjp4GHWsfwYNOL+a929Dt6c/zK+2wyVIoxBx/oEw32JmhwcGiZxI5qoRPm7w45dbg9gemWIvO6dDR3ZDh3S5eTQ+9y63rEN0ukA128VSVWSOmYXP8gTLVYGeGBgf/IpnA/gYP3jzqxqmOwMJxILBVWL41ML06uhsyGHafDrT5O9+FCoE2n4TDTh8sMvDSvjb89CuFEc/D8QcyilRk5r8PcNH5YODbC4Orc/tx1/tOtPgBuwJYJEATgRc+Vw0cVBneDVnn9uPjRh+afQJuf+Bfi1egzR9YbNuuAc0+YENNB7dDIkNiZoyNRc/g5r3biNMuLXCApSwj1yqjyArkWCXkWrt2Q66tdkOSJUwssqDQJod2g5AFoMqdHwgBVQa3QyJDqXP7sWJPCzNjcOzeNLi6dg1+AXT4gVxVhyLLUGQZNgW44Uv2Lt2RwbEJRZUwpiBw6OWOOi9EqMcmcAr0iByF2yGRYYQvT2BmjI1XegZnUyRkqRKEAIJ504WAXYk9szJ6GnWWKsORFdhEN1uVUGiTMb5IhUWRuByBDCN89iUzY2y80jO4smwZzV6BPGtgGnWWIpCtynhlRlHMmZWxdpc4t8ACWQKEJHELMTKcwMnnHahr12BTJGbG4Fj0DCTWIHl1s4bKAgWn3RqKbAquLs/ucSuj4DTq6K2PAHA7JDKc0MnnXj0wCUUDmr3MjJHxFTCQWIPkOoDTbg3jivq+bq67pQfRt/VlX0KidBbs1hyRo6DVK6ALAV2SmBkD40/bQGrdXQfJ7aqEIpsUOsZkxZ4WHG324YRbx8gcBaPy1X4FL3zgX5EkHG3TsKvexwXplFHCJ6GML1Jx0qXBq4OZMTD+pA3EYe86SD62UEVVWWDz28XbmtHi1XCoSYMmBA44/RhX5O9X8KK3XQo/1ZkL1ClTOOwKjrYFCl+WKmNMgQxNCGbGwFj0DGR+pR2/+L/IwyiFLtDm1THv3QY4PTp0BGaiyZIEHcAplw5rvtxr8KK7ZY42+3o9VYEo3cWahMLMGBuLnoHEetcpS8AnTX7Udgi0+wO7xNtVKbB5LhAKe0/Bi9Utc7zVD0e2DLuFpzpT5oqehGJXgH1OPzNjYFynZ3DBKdM2GQiulm3v3BtQAKFd3nsKXqxumaF2GbVunacqUMYLTkJ5pqoIuVaFmTE4Fj2DC4ZueI4CWZKQJQeCqwsBGcCQLAnHW/042uzr9giUWIddZqsKppVZeaoCGUr47zozY0z8aRuc26fhlFuHRwfsKuCTACskKBJQaJVw2qVhRJ6Kei9QW+eNOUAfPtgfpAmBUfkWDsCToTjsCg42+ZgZA+OVnsEdbtbg9Ojo8As0eXSc8QCj8xR8ucQGl1fD6Q6Bz1r8+KzZB58muhyOCfCAWDKPb46wMTMGxys9gzuvUMEplw6PLqAJCTmSQH2HgFXRUecBfDrQ5BXwajpavALji9QuA/Td7TjBbhkyms0nPMyMwfEVMLhsVcHogsCA+36nDx3+wDvTky4tMBtNAnQBSJIEXQiccGm4dJity+PwgFgyg1q3xswYHIueAYSvB4qmibMzyGwy4BYCNllGhx7YAd7vE5CDww6SBL8OdsGQ4XW3HVj0WBwzYzwsehkuej1QNJsshT43zC6j1ScwLEfGFy4dbgD5Vgl5FgmaCJwQPbsii10wZGg9bQcWvVidmTEevlIZoKdNaqPXA0WLHFewYPk0Gzaf8KCmxYedZ3wYapeRrSqhgfZFE3MH81sjSop4MhO9HRgzY2wsemmut01qY60HChdrXGFCiS302BxoJ6PpT2bCd1hhZoyNr1aa6+1daaz1QH3FgXYyov5kpq/bgTEzmY/r9NJcb+9Ko9cDEZldvJnhGjpzYdFLcw670qWghb8rDa4HCm5tRGR28WaG24GZC1/lNNfT0SePbHOGBumDXS7/fSByZ4jw+zDUZAbxZiaIp5qbAy8N0lz0u9JJhUrouKDjLh3b67xYvK055qa3APp0HyIj6U9mgpNfttd5mRmDY9HLAN0dfQJEDtJ3py/3ITKSeDPT0+QXMhYWvQzT2yB9d3hCM5lVXzLT31xR5jFk0Vu9ejUmT54Mh8OBGTNm4KOPPkp1kxKmt0H67vCEZuqJ2TPT31xR5jFc0Vu/fj2WLl2Khx9+GP/85z8xbdo03HLLLfj8889T3bSEmF9ph9AFPm32Yb/Th0+bfRC66HG6NadkU0+YGS5jMBPDFb0XX3wR8+bNw5133onzzjsPzz77LBwOB9asWZPqpiWMLAEQAhIEIMI2v42BU7KpN8wMlzGYidTU1GSYVc1erxfDhg3D73//e8yZMyd0+5IlS7B//35s3ry5y9cU/uHkILaQKD5NC4Yn9fGZGTKa3jJjqCu9hoYGaJqG0tLSiNtLS0tRV1eXolYRpS9mhszGUEWPiIioJ4YqeiUlJVAUBWfOnIm4/cyZMygrK0tRq4jSFzNDZmOoUVqr1YopU6Zg69atEeMTW7duxXXXXRfza5I9ZpIM3R1v8sg2J14+GLmYtjRLRrYKjC+yQhMCVWVW0+wSX11djcrKylQ3I62lY2aS8br1lJnjLh0A8FmzD00eHZIkpSQz0WOlmfi3KVGSmV1DFT0AeOCBB3D//fdj6tSpmD59OtasWYPTp09jwYIFqW5awnR3vInD3nVNkUDg9HSAi20pNrNnJnjMkEcHJEliZgzOcEXvxhtvRGNjI5599lnU1tZi/PjxWLduHcrLy1PdtKSbX2nHL/6vNeI2GcCI3EAx5GJbisXsmQluTm2TAbcIbEXGzBiXoZYsUNcuksuGWZCtKqHFtmZae8Tuzcw02K9bsOuzpsWHnWd8GGqXU5IZdm+exe5N6rfLhmV1GccgorPCuz67G/sj4+CraXBmmbRClAjdjf2RcRhqyQIREVFPWPSIiMg0WPSIiMg0WPSIiMg0WPSIiMg0WPSIiMg0WPSIiMg0uE4vQwUX0da6NTjsXERL1BtmhgAWvYxU5/Zj8bZmePTAPoFH2zTsqvdh5cUFqW4aUVrqKTMsfObC7s0MtLbaHQovgM4d4gXWVrt7+Uoic2JmKIhvcZIoWd0ptW4tFN4gHoFCRlDn9uP3x1Vo9U5mhpKCRS9JYnWn7Kj1YGKRCpcfyFEBQILLL+IOd/gZYEE8AoUyXTAzzhYFeUIPdUH+aEouNp/w4GizDyfcOkbmKBiVrzIz1C8sekkS3Z2i6QKfNPhwpFXDiBwF+5w+SAKYVGKJe3wh/AwwRZJCR6DEOk+PKFMEM9N5fisUSUKrT8Nd/3DCYZdxqEmDJgQOOP0YV+RPWGbIXDimlyTR3Skn2jQISYJPB066NEAAQpJwInRqc9/HF8rsKlZeXICqMisqcmVUlVk5IE8ZL1YX5CmXDrcmcMqtQxcCsiRB77ydmaH+4CueJNHdKR5NAELAKkvw6IAUvF0PnOEb7/gCj0AhowlmJlyHX0eWKkdkRgJCV2zMDMWLV3pJMr/SDpsc6EYBAIsMQAKG5yiwyYAQAgKArbMvh+MLZHbBzHS+D4QmBLItMobZ5VBmAIRyw8xQf7DoJUl0d8rsimycX2yBRZEwPEcBJEASAiNyFY4vEOFsZibna6EuyFdmFCHfqmCYXQ50bQoBGcCwHJmZoX5h92YSRXenBJcw1LVruKBEBSDBrQXerXJ3CKJAZr5b7kdlZVHotpUXK1hb7UZFrg+fu/o3e5MoiL8xg4hjCkTxY24okdi9SUREpsGiR0REpsHuzSTgbu5E8QnPjOJW8dBwPzNDScHfqgTjbu5E8YnOTGurgsXbmpkZSgp2byYYd3Mnik90ZmQJzAwlDYtegnE3d6L4MDM0mNh3kGAOu4JDzb7A3oCagE2RMCxHRlm2JdVNI0pL0ZmRNAkVWRozQ0nBopdg3xxhw8sHXPDrgc1xXX4dLT6BZV+x9fq1nABDZhSdGZ9PwqEmjZmhpOBvR4JtPuHB2AIFp9w6vDpglSUMs8vYfMKDCSXdh5gTYMisojOTDYGKAoWZoaTgmF6C1bo12C0KxhRYML7IgjEFFtgtSq/jE5wAQ2YVnZnybMHMUNKw6CWYw66ETlYI6stu8BzMJ7NiZmgwseglWPSRQn09QaG/wSfKdNGZ0QWYGUoaFr0E6+8Jzf0tlkSZLjozk/M1ZoaShqO9SdCfXeGDwQ8ePcTjhshMwjNTXV3fp997Zob6g78daYRHqBDFh5mheLF7k4iITINFj4iITINFj4iITINFj4iITINFj4iITINFj4iITINFj4iITIPr9PqBx5kQxYeZoXTB37o48TgTovgwM5RO+BsXpbd3pD0dZ8KdIciMmBnKJCx6YfryjpTHmRCdxcxQpsmYiSyvvPIKvvWtb6G8vByFhYU4duxYl/s0NTXhvvvuQ3l5OcrLy3Hfffehqampz8/Rl0MpeZwJZQpmhqirjCl6brcbM2fOxNKlS7u9zz333IO9e/fi9ddfx+uvv469e/fi/vvv7/Nz9OUdKY8zoUzBzBB1lTHdm4sWLQIA7Nq1K+bnDx06hL///e/YsmULpk2bBgBYuXIlrrnmGlRXV6OysrLX53DYFRxtiwxx9DtSHmdCmYKZIerKML91O3fuRG5uLqZPnx66raqqCjk5OdixY0efAjy/0o5d9b5Qd01370h5nAkZATNDZmSYoldXV4eSkhJIYe84JUnCkCFDUFdX1+3XVVdXR3z8AwewsVZFo09CsUVgjsOP5pONaE5ayxMt8o9N9PdnNpn+/fel8PRXOmcm01+3/mF2w/X3++8tMykter/4xS+wYsWKHu/z9ttv49JLL01aG6J/QJUAqiYl7emS78OTER8m849muutrF10mMUNmjPi69QmzG5LM34GUFr2FCxfi1ltv7fE+I0aM6NNjlZWVoaGhAUKI0DtXIQTq6+tRVlY24LYSpQNmhmhgUlr0SkpKUFJSkpDHmjZtGtra2rBz587QGMXOnTvhcrkixiyIMhkzQzQwGTOmV1tbi9raWnz66acAAjPPmpubMXLkSBQVFeG8887DFVdcgcWLF+NXv/oVAGDx4sW46qqrTN1NQObFzBB1lTHr9NasWYPLLrsM9957LwDg1ltvxWWXXYbNmzeH7rN69WpMmjQJN910E2666SZMmjQJv/vd71LVZKKUYmaIupKamppE73cjIiLKfBlzpUdERDRQLHpERGQaLHpERGQaLHpERGQaLHpERGQaLHoGtHr1akyePBkOhwMzZszARx99lOomDYqnn34ahYWFEf/Gjh2b6mZRP1177bVdXs+777471c1KKmY3+dnNmMXp1Dfr16/H0qVL8dxzz6GqqgqrV6/GLbfcgu3bt2PkyJGpbl7SVVZWYtOmTaGPFYUHlWay73znO3jyySdDH2dlZaWwNcnF7A5OdnmlZzAvvvgi5s2bhzvvvBPnnXcenn32WTgcDqxZsybVTRsUqqrC4XCE/g0ZMiTVTaIBsNvtEa9nQUFBqpuUNMzu4GSXRc9AvF4vdu/ejZkzZ0bcPnPmTOzYsSNFrRpcNTU1GDduHCZPnoy7774bNTU1qW4SDcAbb7yB0aNHo6qqCk888QRaW1tT3aSkYHYHL7vs3jSQhoYGaJqG0tLSiNtLS0t7PB/NKC666CK89NJLqKysRH19PZ599lnMmjUL27dvR3FxcaqbR3G65ZZbMHLkSAwdOhQHDx7Ez372M+zbtw8bNmxIddMSjtkdvOyy6JFhXHnllREfX3TRRZgyZQpee+01fP/7309RqyhcPOcB3nXXXaHbJk6ciFGjRuHyyy/H7t27MWXKlOQ2lAbVYGaXRc9ASkpKoCgKzpw5E3H7mTNnTHk+Wm5uLsaNG4cjR46kuinUaSDnAV5wwQVQFAVHjhwxXNFjdiMlM7ssegZitVoxZcoUbN26FXPmzAndvnXrVlx33XWpa1iKdHR0oLq6OqmniFN8BnIe4L59+6BpGhwOR4JblXrMbqRkZpdFz2AeeOAB3H///Zg6dSqmT5+ONWvW4PTp01iwYEGqm5Z0TzzxBK6++mqMGDEiNC7gdrsxd+7cVDeN4nT06FGsW7cOs2bNQnFxMQ4dOoQnnngCkydPRlVVVaqblxTM7uBkl0XPYG688UY0Njbi2WefRW1tLcaPH49169ahvLw81U1Lui+++AL33HMPGhoaMGTIEFx00UV45513TPG9G43FYsH777+P3/72t3C5XBg+fDhmzZqFpUuXGnbtJbM7ONnleXpERGQaXKdHRESmwaJHRESmwaJHRESmwaJHRESmwaJHRESmwaJHRESmwaJHRESmwaJHRIby6quv4oMPPkh1MyhNcXE6ERnCa6+9BrvdDpfLhZEjR+LQoUP40pe+hCuuuCLVTaM0wis9IjKEW2+9FSdPnsSqVavwn//5n7Db7bj88ssj7vPmm29i5MiRaGpqSk0jM9C4cePiOt7n5ZdfxqRJk+DxeJLYqv5j0aN+++KLLzB06FDcc889Ebfv3LkTQ4cO5Rl2NOgkSQr9X1GUiI81TcNTTz2Fu+++G4WFhSloXWzPPPMM3n///VQ3I6ampiacPn0aEydO7PPX3H777fB4PPjDH/6QxJb1H4se9ds555yDO++8Exs2bMBnn30GAKipqcG8efMwbdo0rFy5MsUtJDNZt24dhg0bhoULF+LHP/4xWlpa8Pe//z30+S1btuDQoUNpdWrB0aNH8dRTT+HEiROpbkpMBw4cAABMmDChz1+TnZ2NuXPn4je/+Q10XU9W0/qNRY8G5KGHHoLVasXzzz+P5uZm3HbbbSguLsbatWthsVhS3TwykXnz5uGGG24AELjiu++++yLG8/70pz/hggsuwKhRo1LUwq52794NAAk/FDdRXYsHDx4EgLiu9IDAiREnTpzAP//5z4S0I5FY9GhAHA4H7r77bqxbtw633norGhoasG7durTqPqLMFm83+ne+850uh492dHTg3XffxTe+8Y1un+fjjz/GggULUFlZCYfDgalTp+InP/lJxH22bduG2267DRUVFSgvL8eNN94YKgxBc+bMwRVXXIFPPvkEt912G0aOHImxY8dixYoVEfebNWtW6Krzq1/9KgoLC1FaWgq/39/n5wKAm2++GTNmzMD//u//4rrrrsOwYcPwwAMPxPUYAPDWW29hxowZGDp0KC655BJ88MEHOHDgAMrKyjBkyJC4XocpU6agqKgImzZt6vbnnSqcvUkDdvr0aUyYMAGqqmLTpk2YNm1aqptEBvPoo49i9erV2LlzJ8aMGYOamhpcccUVmDBhAt54441eexW2b9+Oq6++GmvXro15EvnmzZuxYMECjBgxAvPnz0deXh727t2LAwcO4G9/+xsA4C9/+QsWLlyISy+9FNdccw18Ph9Wr14Nt9uNf//736E3emPHjsWQIUPQ0NCAefPmoby8HH/+85+xY8cO/PWvf8VXv/pVAIEis3z5crS2tuJHP/oRAMBut2P27Nl9fi4AGD9+POx2OxoaGnDHHXdgzJgxOOecc+B0Ovv8GKtXr8aSJUtwzTXX4KqrrsLhw4fx6quvYvjw4SgrK8PGjRvjfh3mzJmD+vp6fPjhh/G+3EnFQ2RpwJ577jnoug5N01BUVJTq5pABPfTQQ1i7di2ef/55PPXUU3F3ox8+fBgAYnZtHjlyBPfeey9mzpyJNWvWIDs7O/S5YDfhoUOH8B//8R946KGH8MQTT4Q+P3v2bFxwwQXYuHEj7rrrLtTX16Ourg5erxf/+Mc/Qs935ZVXYtKkSdi1a1eo6F133XV48sknMW3aNNx2222hx+zrcwGA0+nEqVOnkJ2djXfeeQeTJk0KPcYdd9zRp8fYv38/HnvsMXz/+9/HL37xi9B9VVXFf/3Xf+HrX/96v16HUaNG4aOPPurxdUkFdm/SgLz00kt4+eWXsWzZMuTk5GD58uWpbhIZ0EC70Z1OJwDEvP+yZcugqipWrVoVUfAAwGazAQjMsCwtLcV9992HhoaG0L+8vDwUFRWhpqYGALBv3z4AwJIlSyIKbLAg2O320G1tbW04duxYqFAF9fW5AOCTTz4BAPzwhz+MeJx4HuO5556D3W7Ho48+GtGOYHEOn8QSz+tQWFgIr9eLlpaWLp9LJV7pUb9t3rwZTzzxBJYsWYLvfe97qK+vx/PPP49HHnkE5513XqqbRwbzgx/8AKtWrcKuXbuwadOmfk1IESJyNMfj8eCvf/0r5s6d220B9Xg82LJlC1wuF8aOHRvzPnl5eQACV01A4IoqXPBK89xzzw3dtn//fgghIopVPM8V/nw33XRTvx7D6/Viy5YtmDdvXsTjAgiNLUYX5b6+DsGfdfiykXTAokf9snv3btx77724/vrrQ+MRixYtwm9/+1ssW7YsbdfoUOYaSDd6cXExgMC6s4qKitDtNTU1cLlcOP/887v92mPHjsHlcmHp0qWoqqqKeZ9x48YBCFzplZSURDwHcPaKLLyABK8Kw2+L57mCjzFkyBBUVlb26zGOHj0Kl8sVc/borl27IMtylzewfX0dmpqaYLVauxTTVGPRo7idPHkSc+fOxYQJE7Bq1arQO7ni4mLcfffd+PWvf41HHnkkrrU9RD0J70Z/6qmnsHz5cqxevbrPXx/8w33s2LGIAtfe3g6g56uR1tZWAEBlZWXE+FYs+/fv73JlBASK07Bhw0LFN3hbSUkJhg4d2q/nCj5G9PPF8xjB7z+a1+vFn/70J4wePTqiyzee1+HYsWPdXmmmEsf0KG7Dhw/HgQMH8M477yArKyvicz//+c/hdDpZ8ChhorvR77vvPqxfvx6HDh3q82NMmTIFWVlZoXVxQaNHj4YkSTF3RPH5fACA8vJySJKEN998s8t9dF1HY2Nj6P8HDx6MWfQ++eSTLmvdPv/8cwwfPjzitr4+FxDoPoz1fPE8xogRIwCgy3q6F154ITQrOyje12HPnj2YPn16zM+lEq/0iChtJaob3Wq14vLLL8fWrVvx5JNPhm7Pz8/H7bffjj/+8Y/w+Xy47LLL4HK58NFHH+Giiy7Cww8/jNLSUtxyyy1Yt24dbrrpJlx11VUQQuDYsWPYtGkTnnnmGVx99dWoqamB2+3uUoQ0TcPBgwdx//33R9xeUVGB9957DytXrsQ555yDc889F1OnTu3TcwFnuyajn6+v7QWAIUOG4Morr8S6deuQlZWFCy+8EO+//36oCAaLXryvw+7du+F0OnHttdf26fUZTCx6RJSWEt2Nfvvtt2Pu3LmoqamJmHzxzDPP4JxzzsEbb7yBLVu2ID8/H1OnTsU111wTus9vfvMbTJw4EX/+85/x05/+FDabDRUVFbj11ltx8cUXA4g9RgcAn376KTo6Orpc6T388MM4cuQInn/+ebS2tuKnP/0ppk6d2qfnCn++WLul9PUxgECX5YMPPogNGzbgzTffxKWXXornnnsOd955JyZMmNCv12HDhg0YMWJEn7poBxsXpxORKei6josvvhhXX301fvazn6W6OYbV0dGByZMn48EHH8SiRYtS3ZwuOKZHRKYgyzIef/xxrFmzhkcLJdEf//hHWK1WfPe73011U2LilR4REZkGr/SIiMg0WPSIiMg0WPSIiMg0WPSIiMg0WPSIiMg0WPSIiMg0WPSIiMg0WPSIiMg0WPSIiMg0/j/soyt/qZPavQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x     = np.random.normal(3,3,100) # generate 1000 points\n",
    "xbar  = np.mean(x)\n",
    "xstar = x-xbar\n",
    "\n",
    "epsilons = np.random.normal(0,2,100)\n",
    "\n",
    "y= 2+3*x + epsilons\n",
    "\n",
    "\n",
    "fig,axs =plt.subplots(1,2) \n",
    "ax=axs[0]\n",
    "ax.scatter(x,y,alpha=0.7)\n",
    "ax.set(xlabel=r\"$x$\",ylabel=\"y\")\n",
    "\n",
    "ax.axvline(0)\n",
    "ax.axhline(0)\n",
    "\n",
    "\n",
    "ax=axs[1]\n",
    "ax.scatter(xstar,y,alpha=0.7)\n",
    "ax.set(xlabel=r\"$x^{*} (centered x)$\",ylabel=\"y\")\n",
    "\n",
    "\n",
    "ax.axvline(0)\n",
    "ax.axhline(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By centering, all of the (x,y) pairs in this particular case move to the left.\n",
    "\n",
    "\n",
    "\n",
    "In our original SLR, we can compute the MLEs for $\\beta_{0}$ and $\\beta_{1}$ as \n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{\\beta_{0}} &= \\bar{y} - \\beta_{1}\\bar{x}\\\\\n",
    "    \\hat{\\beta_{1}} &= \\frac{ \\sum_{i=1}^{N}(x_{i}-\\bar{x})(y_{i}-\\bar{y})}{\\sum_{i=1}^{N}(x_{i}-\\bar{x})^{2}} = \\frac{ \\sum_{i=1}^{N}(x_{i}-\\bar{x})y_{i}}{\\sum_{i=1}^{N}(x_{i}-\\bar{x})^{2}}\n",
    "\\end{align}\n",
    "\n",
    "Our new SLR model where we center $x$ simplifies the above estimates\n",
    "\\begin{align}\n",
    "    \\hat{\\beta_{0}^{*}} &= \\bar{y}\\\\\n",
    "    \\hat{\\beta_{1}^{*}} &= \\frac{ \\sum_{i=1}^{N}x^{*}_{i}y_{i}}{\\sum_{i=1}^{N}(x^{*}_{i})^{2}}\n",
    "\\end{align}\n",
    "\n",
    "Above we used the (very handy) fact that \n",
    "\n",
    "\\begin{align}\n",
    "   \\bar{x^{*}} &= \\frac{\\sum_{i=1}^{N}(x_{i} - \\bar{x})}{N}\\\\\n",
    "               &= \\frac{1}{N} \\left( \\sum_{i=1}^{N} x_{i} - \\sum_{i=1}^{N}\\bar{x}\\right)\\\\\n",
    "               &= \\frac{1}{N} \\left( \\sum_{i=1}^{N} x_{i} - N\\bar{x}\\right)\\\\\n",
    "               &= \\frac{1}{N} \\left( N\\bar{x} - N\\bar{x}\\right)\\\\\n",
    "               &=0\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Are $\\beta_{1}$ and $\\beta_{1}^{*}$ equivalent?\n",
    "\n",
    "Shifting all of our collected xvalues by the same constant does not change the relationship between Y and X.\n",
    "We can see that this is the case two ways. \n",
    "\n",
    "###### Algebraic\n",
    "\\begin{align}\n",
    "    Y_{i} &=  \\beta^{*}_{0}+\\beta^{*}_{1}x^{*}_{i} + \\epsilon_{i}\\\\\n",
    "          &=  \\beta^{*}_{0}+\\beta^{*}_{1}(x_{i}-\\bar{x}) + \\epsilon_{i}\\\\\n",
    "          &=  \\beta^{*}_{0}-\\beta^{*}_{1}\\bar{x} + \\beta^{*}_{1}x_{i} + \\epsilon_{i}\\\\\n",
    "    \\epsilon_{i} &\\sim \\mathcal{N}\\left(0,\\sigma^{2}\\right)\\\\\n",
    "\\end{align}\n",
    "\n",
    "We see that our $\\beta_{1}^{*}$ is attached to $x$ in the same way as our original SLR.\n",
    "This then means $\\beta_{1}^{*} = \\beta_{1}$\n",
    "\n",
    "###### Visually\n",
    "We can plot the points $(x,y)$ and $(x^{*},y)$ and see that the estimated **slope** remains unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x1c917fbe0>"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAEjCAYAAACmQt93AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABNzUlEQVR4nO3dd3wT9f8H8NflMtp008koIJICBRlWoaCIgCyRpSCCrLLLLFAUFBAEKYVSQMCClCGCX6lsEfGngijKUhnKLEIrs6F7pG2Sy/3+iA29Nt3jMt7Px8PHw37uc8knTV+8k/vcfY5JT0/nQQghhNgBidgDIIQQQmoLFT1CCCF2g4oeIYQQu0FFjxBCiN2gokcIIcRuUNEjhBBiN6joEUIIsRtU9AghhNgNqdgDIOXjvv2+4Of0kPoijaR2xMfHQ6VSiT2MGmUPr7GqbPl3RJkWB33TI4QQYjeo6BFCCLEbVPQIIYTYDSp6hBBC7AYVPUIIIXaDih4hhBC7QUWPEJF4f/UVmORksYdBiNWQb91a5cxQ0SNEBLKvvkKjlSvh/NJLYM+cEXs4hFg82VdfwXHOnCpnhooeIbVMcuMGHMPCjP//4AGcXnsNsl27xB0UIRZMcvOmMDN9+0L2+eeVeixakYWQWqDW6LEzXoO01CwsnD0SLjk5TzayLLjWrcUbHCEW6ElmsrFw9ghhZqTSSmeGih4hNUyt0WPW6Qzkcwa8s3EBfO/cFGzPXbkSBip6hJgIMvOJmcxERsLQpk2lHpsObxJSw3bGa5Bv4PHqib3ocfKAYNvFHoMQpuqHqEuZUGv0Io2QEMtSkJk+J/ahx0/7Bdsu9RiEsID+lc4MFT1CaliShkPTxOuYsnWxoP2mXxMsHLkY/2p4nFFrMet0BhU+QvAkM1O3fiBoj6/7NBZUMTN0eJOQalQwD5Gk4eCrZDFKpURDXoOQqKlQaPNN/TRyR4ROiIZM6QwAYBkG+QYeO+M1CG/jKtbwCal1JWVmzOppxTIzfco6oIqZoaJHSDUxzUMYeLAMgzvZHC481mLHtvnwepQo6PveyCXIatwEdQq1sQwDdS5Xu4MmRESlZuZhgqDv+yMWI0/VDA6F2iqTGTq8SUg1KZiHYBkGgDGQvY/uhNfRw4J+5/sNB/PWEMiKpI/jefg4srU1XEJEV+7MvDYMGPYmZCwjaK9MZqjoEVJNkjScKbwAEBB/CZN2Rgj6cM88g4At0ZjS0gVyCQ+O543tPA+FhMEolbJWx0yImIpl5tbl4plp1cqUGYWEqXJm6PAmIdXEV8niTrYxxM5Z6Xg/ejpkep1pe56TCx5t2goPBwf4AJj3tBanOA+oczn4OBrnMnyUFEliP4plZvU0YWaUzsbMODrCB8Cajm7YGa+pUmas5pveli1b0KlTJ/j7+8Pf3x89evTAd999Z9rO8zwiIiLQvHlz+Pn5oW/fvrh27ZqIIyb2ZpRKCYWEgYHjEL5hLnwf3xdsj5wcgWmPPE1nm3kqgPA2rlgZ7IHwNq7VXvAoM8TSlZWZlZMjMC3Jy5QZH6W0ypmxmqJXr149LFmyBCdPnsSJEyfw0ksv4e2338bff/8NAFi3bh02btyIyMhIHD9+HN7e3hg0aBCysrJEHjmxFz5KKdZ0dMOc49sR/MdxwbYDfcfgdMc+prPNagNlhli6gszMPrGjWGYOvjoav3V6tdozYzVFr2/fvujRoweaNGmCpk2bYuHChXB2dsb58+fB8zxiYmIQFhaGAQMGIDAwEDExMcjOzsbevXvFHjqxI3UvnEXP2JWCtmuqttg64l0AtXuGJmWGWIO6F8+h15ZIQdv1pm0QO3IegOrPjNUUvcI4jsO+ffuQk5OD9u3bIzExEUlJSejWrZupj6OjIzp16oSzZ8+KOFJiT5jHj6EcNw4M9ySgmc7uWD57PfQyOQDxztCkzBBLZDYzTm5YPqfmMmNVs+ZXrlxBz549kZeXBycnJ+zatQstW7Y0hdTb21vQ39vbGw8fPiz1MePj42tsvNVLeIaS9Yy78qzqNXIcAqZPh+TRI0FzxISluO3gBkl2Dgw8IJfweJFNQ3x8EoCyX6NKparSsGwhM1b1d1Ahdp5pjkPAjBmQFPl7i5i4DP84uJeYmbKUlRmrKnoqlQq//PILMjMzcejQIYSGhuLIkSNVfkyrcEo4wWs1466k+Ph4q3qNiogIOJw/L2jLmz0bY8OHQVrC2Wa18RqtPTPW9ndQIXaeacWKFXA4d07QJ2/WLIydW3JmqoNVFT25XI4mTZoAANq2bYs///wTn3zyCcLDwwEAjx8/hr+/v6n/48eP4ePjI8pYif2QHj8OxUrhPJ7+hReQ/9578JFKRV1WjDJDLJH0xAkoIoXzePpOnZD//vs1nhmrnNMrYDAYoNVq0ahRI/j6+uLEiROmbXl5eTh9+jQ6dOgg4giJrWMePIDjhAlg/rtgFgAM3t7QbN0KSC3vMyVlhohN7MxYXipLsHjxYvTs2RP169c3nWF26tQpxMXFgWEYhIaGIjo6GiqVCk2bNkVUVBScnJwwePBgsYdObJVOB+W4cZCkpJiaeIaBJjYWvJ+fiAMzoswQi6PXGzOTnGxqMmWmbt1aGYLVFL2kpCRMnDgRarUarq6uaNmyJfbu3Yvu3bsDAGbOnInc3FzMnTsX6enpCAoKwv79++Hi4iLyyImtcli6FNLTpwVt+fPng+vSRaQRCVFmiKUxm5l582o1M0x6ejpfdjciNvftwknv9JD6Io2kdlj6CQzSo0fhNHy4oE3XvTs0X30FSMo3a2Dpr9ES2PLvyN4yrd66Fao5cwRtum7doNm7t9yZqQ5WPadHiBiYhAQoQ0MFbYZ69ZC7eXOthpcQa8EkJuKpJUsEbYZ69ZD76ae1nhlKKCEVkZ8PZUgImIwMUxMvlUKzfTt4Ly8RB0aIhfovM9LMTFMTz7LQbNsmSmao6BFSAQ4LFkB64YKgLe+DD8DRGY+EmOWwcCGkf/4paMtbvBhccLAo46GiR0g5yfbvh2LLFkGb7tVXoZ02TaQREWLZpAcPQvHpp4I2sTNDRY+QcpDcugXHGTMEbal1/fFv1McAw5SwFyH2S3LrFpTTpwva0vzEzwwVPULKotFAOWoUmOxsU5NWKseCmesx8xpjutcXIeQ/ublQjh4NptBtqrRSGRaEfSx6ZqjoEVICtUaPqEuZ+GtMGNirVwXbNocswO2mz9Tq/fEIsXQFmbk8ZibYK1cE29YPm4t/mrYWPTNU9AgxQ63RY9bpDLju+R9e/L84wbafXngN3/Q0XqNXm/fHI8SSFWTGJe5LdP6ueGYOdhsKQPzMUNEjxIyd8RrUvXMd02M/ELT/4/cU1k3+yDQnIdb98QixNDvjNaibcAMztiwStN/2bWxRmbGaZcgIqS5qjR474zVI0nBwkgIAgxw9D1/lk9uYpCdnYFH0dDho80z75ckdMHXSWnAOTmBhDK9CwmCUSlnSUxFiE8qVmZRMLFo9TZCZfJkCUyavA+fgBORoLCIzVPSIXSk4BJNv4KHjeFxJ04HhgVaeMtzJ5nAhWYc1wa6YsPF9+D+4Ldh3/fglaPFiazjL2Rq71xchlqZKmZnwJDM3H2YjwEcuemYorcSu7IzXIN/Ag2UYJOToAd64yvu9bA5N3GTIN/C4FvUp+h8/LNjv266DcbL7YKxp6UJFjtiVcmVm9afo/+MhwX7Hug7GT92HmDITr0yCSiXevSULUHqJXUnScGD/m1vINwCM6f+N6643u/03+mwQrhF4r3EznJ/3EdY840YFj9idMjNz5wr6rC+emXMWmhk6kYXYFV8lC+6/m1cqJADP8+ABKCQMnLMz8P7q6ZDptKb+vLMz3PbuxswOvhYXXkJqQ2mZccrJxPtR04pn5qtdFpsZKnrEroxSKaGQMOB4HvWdWIABGJ5HAycJZm18B3XVdwX9cz/+GIamTUUaLSHiKzUzG0rIjAXfDsryyjAhNchHKcWajm7YGa+BOpdDO08pAAad4j7FC+d/EPTNnzAButdfF2eghFiIkjOzBS+e/17QN3/8eIvPDBU9Ynd8lFKEt3kyoc6eOQOnbZGCPvp27ZC3bFltD40Qi1QsM2fPwml7kcy0bYu8jz6q7aFVGB3eJHaNSU6GcuxYMNyTFSJ4Nzdotm8HFAoRR0aIZWJSUoyZ0T9ZP5N3c4Nmxw6ryAwVPWK/DAY4TpwIyYMHgmZNTAz4xo3FGRMhlqwgM/fvC5qtKTN0eJPYlcIrS7y5dyNeOX5csD1/xgzoX31VpNERYnkEmdn3CV758UfB9vzp060qM1T0iM0qHFZfJYtXGyjw0cVs5Bt4PPv3aXTdFi3or+/YEXkLF4o0WkLEV1pm2l05g65bVwv664ODkbdoUQmPZpmo6BGbVHjpJJZhcCebQ9wtDfycJKiflYJ5a8PA8gZTf4OnJzRbtwIymWn/wuEXe+kkQmqavWTG8kZESDUovHQSYLydiYbjoc7UYv3amfDISDH1NTAMcmNjwderB8B8+C8k67Cmo+WtLkFIdSk9M2Gok55s6mtgGORu2QK+fn0A1pUZOpGF2KTCSycVcGAZTD2wDq2vnhO0/zR6JvRdu5p+Nhd+sW98SUhNKykzUw58jNZXzwraT46aAX23bqafrSkzVPSITSq8dFKBQTd+xrRvtwjaLrR+Af5L5wvazIVf7BtfElLTzGVm4M1fMP3bTwVtF5/phAbL3hO0WVNmqOgRm1R46SQAqKO+jwWfvCPok+npA7ddsfBxEV5bZC78Yt/4kpCaZi4zCzfOFfTJ9PSBq5VnhooesUkFSycF+8jRRKFH5IaZcM1ON23nWRaSndvh2bBusX2Lht8SbnxJSE0rmpkVG8KEmZFIjJlpVK/YvtaUGcuaYSSkGhUsneQwbzkUVy8ItqW8+z5kL7xQ4n6F1xqkm8USe2HKzPwIKK7+KdhmK5mxvBERUknmTpmu9/03UGzaJOh3OqgrlrcbiWiNvsRQFl1rkBBbZDYzPxyFIiZG0O/Ms13x0bOjbCIzVPSITTB3yvSjv+Kxaf40Qb8k7/pYPS0KeWCwM15jFSElpCaYy8zDv29h87ypgn5JXvWwevoqm8kMFT1itQp/Sr2ZoYeB56GUGSfOHbX5eDdqGtisLFN/nVSGj2avR7aLO1jAIs8sI6QmlZYZB53WbGaWz16PLBcPm8kMFT1ila6m5GPMT2nQcDwcWAY5Wg5aA+Ai58CBwcpdS9A04Zpgny2j5uOmqg0Ayz2zjJCaUlZmInctgerOVcE+sSPn4UZAWwC2kxkqesTqqDV6jDmZhpR8AyQMgzy9ATk6QMcDGs6AMX8ewVsn4wT7/NLxVRzoPRIsLPvMMkJqQlmZGX3hGww7uUewz6/BvbG/zyibywwVPWJ1dsZrkKszhhcAGIYBwMPAA80e3sKqLxYL+nNNmqDejo0IfsRa/JllhNSEsjITtfsDQX+uSRP47fgEwUlSm8uM9b8CYneSNBwcpBLkankUrAHBM4CbToMvt8+CkzbX1FcrVyD/s8/g5euBcF9xxkuI2ErLzP92zBZkRieTI2/HDnj51UG4nzjjrUl0cTqxOr5KFnWdJJAAKFgDguF5bIz7EM0e3Rb0/SZsKQzPPFPrYyTEkpSUmQ1fLUXzh/8I+n4TthSG1q1rfYy1hb7pEaszSqXEhWQdmnsAD3MMyOMMGPf7Pgz7/WtBvx9fHgS/iaMQdSnT4m93QkhNMpeZsX/sx/DzhwX9jncZCN9Jo206M7bzSojdKLr6Q5u71zA5brmgT9JTAXBevxpLL+VYxe1OCKlJRTPT+t51hO75SNBH3VgFpw3RNp8Zqzm8GR0dja5du8Lf3x9PP/00hg4diqtXhafX8jyPiIgING/eHH5+fujbty+uXbtWwiMSa1aw+sPKFhJMWj4Vkvx80zbeyQnKLz/HkVSp1dzupCZQZkhhhTMzuWhmlEo47tllF5mxmqJ36tQpjBs3Dt999x0OHz4MqVSKgQMHIi0tzdRn3bp12LhxIyIjI3H8+HF4e3tj0KBByCp0sSWxITwP5bRpYO/cETTnrl0LQ7NmVnW7k5pAmSHF8DyU06eDvS2c+7anzFhN0du/fz9GjBiBwMBAtGzZEps3b0ZycjLOnDkDwPiJNSYmBmFhYRgwYAACAwMRExOD7Oxs7N27V+TRk5og37QJsq+F83j5Y8dCN2QIAOu63UlNoMyQouSbN0N2WDiPlx8SAt2bbwKwj8xYTdErKjs7GwaDAe7u7gCAxMREJCUloVuhu/k6OjqiU6dOOHv2bAmPQqxV5i9nIF+wUND2QNUKdxd8aPrZmm53UhsoM/Yt49RZyIpk5mHTlri7cKnpZ3vIjNXOTM6bNw/PPPMM2rdvDwBISkoCAHh7ewv6eXt74+HDhyU+Tnx8fM0NsloJ/+isZ9xPpOQDB5OkSNEx8JTxGOirh6ei5P4lvcZMdTpajh4DltOb2jIcXTAlJBIZv6Rg3tMPTY873df4nKk6BnX+e86M+6nIqM4XVgVlvY8qlaranstaM2ONf+vlU3amK5qZkmSq0xE4egyket2TNgdnhI5dWauZqY33sqzMWGXRe++993DmzBkcO3YMLFu1r93V+Y9KjTp1X/Cj1Yz7P2qNHosLVnSXMkgz8FifxJR4Vlh8fLz512gw4NH0wfBNfSRonjHyI1zxfRqNlVKc4jwQrjKuBK8CENyqJl5R1ZX4GmuAtWamNn9Hta6MTFc0MyX6LzN+KcIPMrWdGUt5L63u8Ob8+fOxb98+HD58GI0bNza1+/oal9t4/PixoP/jx4/h4+NTm0MkZuyM11TLWWGKtWvR7MxxQdv6rqNwpO0r0Bpsb9K9OlBmrFO1ZWbdumKZ2fDySHzdroddZsaqit67775rCm9AQIBgW6NGjeDr64sTJ06Y2vLy8nD69Gl06NChtodKiqjoWWEp+UDUpUzMPZ2GqEuZUGv0YE+dgmLZMkG/043b4t1+s5Gl5SEBb3OT7lVFmbFeFc2MWqM3n5mlSwX9zjRuY9eZsZrDm+Hh4dizZw927doFd3d303yEk5MTnJ2dwTAMQkNDER0dDZVKhaZNmyIqKgpOTk4YPHiwyKMnvkoWd7KNIdZyBtzL5pCnNyBPL4O6yN2Y1Ro9Vvwjh1SpNV0ge+fWQ+yYPxaMwWDql+LkjrdHRyGPkYHhgdR8A3gDb1OT7lVBmbFuFc1M0RvC3rn1EDveGyfITKqTG94eHYVcidxuM2M1RS82NhYAMGDAAEH7u+++i/nz5wMAZs6cidzcXMydOxfp6ekICgrC/v374eLiUuvjJUIFyyBl6ThcT+PA8cZwGnges05nmOYp1Bo9JvychhvZErjo9WjgzMIBPOasCYNUrRY85pxxkXjoWRes4ckhCwlT/LntFWXGulU0M7czdHCQSoSZ+e+DToHZ41bigWc9u84Mk56ezpfdjYjNfbtw0js9pL5II6m8gnD+k6mHA8ugvhMLB6kEHM8j2EeOUSolZp3OwKUULXK0HFiWhQTAqh82YNS+jYLH2v36FCzuOwPp+QYwDAMegIecQSNXKYJ95Ahv4yrKa6wIS5nYt2S2/DsqT6Yrkpl8zriYtATAyh82YvS+DYLH+mJQKD54baZombGU99JqvukR6+ejlCLATQoHqXAquWCeYme8Blk6Dho9jxw9IOMN6H3jV4zY/4mgf0aHTtjz1kzkpT8JrwRAA2fW7ibliW2rUGZ0PKQSoNeN3zByv/BDYkaHTvhyWBhlBlZ2IguxfqWt+JCQqcf1NA48z8MAwDv5ITZtfweSQv1T3bwwdUI05ge542lXKRxYHh5yBi3rSCFnJXY3KU9sX3kzwwPwSnlkNjPTJqymzPyHih6pFQVnlt3J0OHfLD1y9cZPloVXfLibY5y3kEokcOV1+GJnOLxy0k2PwUkkWDFrHZLcvHD0Xj62vOSB1p4KNHJ9El5bWz2C2K+KZsadMWbGO+fJ2qrGzKzFIzdvysx/6PAmqXFFzyzzdZTgUY4B7X2kaOwqM92vq4FSgmtpxon6JV+vRcc7FwWPs3PoLFxuFQyd3oBj/+YiScNB5SoBwEDDGT+t2tq9v4h9qkxmPjj8MTrdviB4nM+HhuFyq46UmULs55US0RS9yFYpY9HQVQIvB+PPqy5lwVfJwttRgmbuLAJ//R7Tjn8meIzzbV9C3KDJyNMbcCVNBzcZg39zDKZPqrZ0vy9CKpOZqT9sEzzG7207Y8+gUMpMEXR4k9Q4cxfZ6jgeBxPycUatxb85BpxRa3ElTY9Gj+8iats8QV+1px9WTI8CL5HgXg4HhgcauhjDaov3+yKkIplpmHwPq7bPF/RVe/ohYvpqyowZ9lfmSa1zkgK3MnTQGgCFBKjvxOJeDgeWefJJVsfxUKdrsHb1DDhrntzLjZdKkb99OwLd/KDO5ZCSK0FDZxZy9snnNXs7+4zYvvJnJhdrVk+HS06maV9eKkX+tm0IdKfMmENFj9QotUaPK2l6ZOYbwDMMcnke6VoDHCSAyl0GLWfAnQwdkvJ4rImLQGDiFcH+eUuWIO+554F4DXgecJIx0Bt4yAudbGZvZ58R21aRzER/FYGWCUUys3gx8p5vT5kpARU9UqN2xmvAMwxaecpwL5tDvgGQMYCPo/E07KupeqRreQz+4ygmnfpSsO/Vzr0hGzNRMKHPg8fNDA4BbsZ5Dns8+4zYtvJm5o0/v8XkX/4n2Pfqi70gC5lEmSkFFT1SowrmJliWQRO3J4dXvBXA6SQd0rU8mj66jS17Fgn2S/RqgHFvfQjFL+kw8DyUMuOnUkcpi2buAAMGjZwldnn2GbFt5cnM00l3imXmrlcDjBu2lDJTBvt81aTWFF40twDH82jsKsflFC1c9bmI+2w2nLW5pu15UjmGjo7Gdd4Z/GMt5FIGLT0Y06oUjlIWjZwlWBnsUeuvh5CaVp7MfLVjFlzyn5yIks/K8NaYaFyjzJSJzt4kNWqUSgmFhDGtKFH40IqcZbFh30d45qHwbsozB83Hn/6BYCWAHkBGPo9/s5/cJd2e5yOI7SsrM+v3LTebmfP+LSkz5UBFj9QoH6UUazq6IdhHDm8FAx3Hw11mnLeY9sd+DDt9QNB/97N9sa3TEChZ4xlmDiwDMEBKnvH2KPY+H0FsX2mZmfrHAQw/vV/Q/4tnX8XWF96Ei4yhzJQDFT1S43yUUoxSKZGm5SFjGSRrgaTzl/HWxg8E/W75NcHc4R/AVc6YbnfCMgycpYCz1DgfEewjt9uLaon9KCkzwzYK5/Fu+TVB+PDFcJUzgjusU2ZKRr8FUisKrzCh1GRh0eppUOjyTdu1Cgd8u2IzXvPxwoHEPOj1xkM7PM+DlTB482lHLH6e5iOI/SicGcfcbCyKnl48MxGb8ZqvMTM8zxvvoECZKRV90yO1wrTCBM9j5qb30OBhgmD7R2OWYCfbEACPQDcWLiwPRykDN4UEberIMKUl3dSU2BdBZmLeQ4MHdwTbl49ejJ3SJ5lxU0goM+VA3/RIrSg4I23gsV3o8ttRwbbtHd/Ahjb9IMnikJDNobkri25eeihcne3+9GpivwoyM+C73Xj5t28E23YEv4H1bfsLMtOvoYPdLiJdEfRbIbVilEqJ7DN/YOKOjwTtF+s1w9SB70GXz8NVwUDCAw/zeLRyAFbQ6dXEjo1SKZF15g9M3LFM0H6pXgCmDCqeGWe5BItr+O7ntoCKHqkWao0eO+M1SNJw8FUW/6Tpo83GyvUzION0prZMhROGjlmDXLkDGAAaHQ8XuQQ6A5CqY8w8CyG2S63RF8vMqg0zIdMXycxo85mx17U0K4qKHqkytUaPqafScDNDj/R8HjzP43/xGnze1QOBngqA56GcPBmye3cF+015eylueTcCA4AHoOONJ67IJAzqyHizz0WILVBr9MXa+hxNFmYmNBSyu/8K+kwd/iHifRqbzQxdh1c+dCILqbJPrmTjz2QtHmoM0HA8cg1AYjaHESdSodboIV+/HrJjx4T7vPQ29rftBfa/L3QMYEwxAzR2ZjHQt/g/CoTYik+uZBdrE2RmwwbIvv1WsD3mpeHY1653iZmh6/DKh77pEbPKOlxZ2KlHWuTqAZ4HuEJf0NS5Bvy07yeMX7JE0P/iU89gbv9w/HdVAhQSQGcAvB0ZvPW0I6a0dEHG/fQaemWE1IyKZqYoA/9fZvafxPjFiwXbLjZ+BnP7zy01M3TiSvnQb4kUo9boBau038nmcCFZV8oFrjz0BpgCWVD3lGkpeH3tVDDck7mGdKUbJk6IBi+XgzUAkACOLNDZW46NL3qYHj+jRl8hIdWrMpkpygDAMT0Vr6+ZUiwzkyZEw1BGZkj50OFNUkzhi2KBsu+0/KKfHAYYY1wQZYmBw65d78IzNUnQd8roCCR4NIADC7AsUEfB4GlXKYWXWLXKZKYoiYHD5yVk5k4dykx1oaJHijFdFFtIaXdantLSBV4KYdv7329Gj5unBW2be4/H9627AgCkEgkcWQaOUgl6N3Q0hVet0SPqUiZW/iND1KVMsxP+hFiaymSmqPd++BS9bvwmaNvce1y5MzP3dBplphyo6JFifJWsaYX3AqWt0u6jlOKtpk5wlgJSBuh58zQWf/eJoM/1ls8jatBM5OgMyNEZYOB5gGGgN8A0AV9wiOiMWouHeRKcUWsx63QGhZhYvMpkpqglxzYKfjZmJqzcmfk3x0CZKQcqeqSY0m5tUpI3mzhCwTKon6HG55+/A0mhfwBSXOogJCQKaRwLR6lxBXiNnoeTFBjU2MH0D0BFDxERYikqk5miimZm7JhVlJkaQAeESTEFtzbZGa+BOpcrc1kjtUaPjy5mI8DJgMiocPhkp5q2GRgGM8atxENXb7AcAwPPw0kqgYHnoWAZTGnpbOpb0UNEhFiKymSmJAWZeeDmQ5mpAVT0iFk+SinCy7mkUcGnzZkH1qPz7T8E27a+Pg2nAl8AGAYtPFjcz+GgNQByCYP23jLBPwol3TGaLrol1qCimSnJ1kGUmZpERY9UWZKGQ8c/TmDowc2C9h+adcKibpOg4A2QsiwcpBI87WY8os7xPBq7ygT9R6mUuJCsQ76h8oeICLEGSRrz38R+DOiIRd0rlxmWYSgz5UBzeqTKmmU8QPj6cEHbPTdfjB8dCR0jgZSVoLEzW+Z8R+E7RtdzMNDNL4nNCsh8WKztvpsPxo+pfGbohrHlQ78ZUjVaLaaumAaHnCeXk+slLEJCVkHv6YU6UgYdfORYFORarvmOgkNE8cokqFS0YjyxQVotpq2YhrmjPhc0h4yJgs7Tu9KZIeVToaL3448/olu3bmAYWgGfGDksWADFxQuCttVvzAbXsSNasZL/DslI7TaYlBlSlMPChVBc+BMYJWzXd6LM1IYKHd4cPHgwAgMDsXDhQvz11181NSZiJaQHD0Lx6aeCtrPPdccPQyZC/l947X1+gTJDCpMeOgTF5s1mt1FmakeFit7u3bvRoUMHxMbGokuXLujUqRPWr1+Phw+LH58mtk1y6xaU06cL2gwNG8Jr52YE+ypofuE/lBlSQPLPP1BOm1bidspM7WDS09MrfOOyrKwsHDx4EF999RV+/fVXMAyDzp0746233kK/fv2gVNKnlOrmvv2+4Of0kPqijEOt0eOLv1MxdlJ/+N+5bmrn5XLkfPcduHbtquV54uPjoVKpquWxLIG5zAQFBWHs2LGUmVLYwt+BKTOTB8D/9jVTO7PmiqCfWJmuLZbyXlbq7E0XFxeMHDkShw8fxl9//YVFixbh8ePHCA0NRUBAACZNmoSTJ09W91iJyAqWPHo+coGg4AFA3kcfVVvBs0XmMpOamkqZsXEFmXkucqGg4BHxVPmSBY7joNPpoNVqwfM8HBwccPLkSQwcOBCdO3fG1atXq2OcxALsjNeg84n9ePXHOEH7X11fg3b8eJFGZX0KMqPT6SgzNq4gM31/3CNo//vlviKNiFSq6GVkZOCzzz7Dq6++irZt22LlypVo1qwZdu3ahevXr+Pq1av4/PPPkZGRgalTp1bbYH/99Ve89dZbaNGiBdzd3bF7927Bdp7nERERgebNm8PPzw99+/bFtWv06aq6sFevYuanCwVt9+o9hS3TI4AKnJ1oj6vCm8vMU089RZmxcZJr1zBzyyJB2726jY2ZqQB7zExNqdBs6ZEjRxAXF4fvv/8eeXl5ePbZZ7FixQoMHjwYHh4egr6vvfYaUlNTMWfOnGobbE5ODgIDAzFs2DBMnjy52PZ169Zh48aN2LhxI1QqFVauXIlBgwbh/PnzcHEpfisPUgHZ2ZgaMQ0O2jxTU75cgQ9nr4evp1u5H6biN9u0bqVlJjk5WTDHQZmxMdnZmBoxFQ75uaamfLkCH87ZAF8vd+DfrHI9jL1lpqZV6Dc2cuRI1KtXD5MnT8awYcMQEBBQav+WLVtiyJAhVRpgYT179kTPnj0BAFOmTBFs43keMTExCAsLw4ABAwAAMTExUKlU2Lt3L0JCQqptHLZKrdFjZ7wGSRoOvsonF8Oqc3TQjJmGNv/eEvRfP34JHjzVAnMrcHp1aavC2+I1SaVlJjk5uVh/yox1KTMziUUyM26xKTPL/ixf0bO3zNS0ChW9AwcOoEuXLuW+0DYoKAhBQUGVGlhFJSYmIikpCd26dTO1OTo6olOnTjh79iwF2IzCgXWSMvgrVQdGwgg+Tb7f1hl/rPwU0344KNj35+5vIGPocKwpZSV5c+xtVXjKjG0pb2Z+X7kF04tk5pdurxszE+BEmRFRhYreyy+/XEPDqLqkpCQAgLe3t6Dd29ubrokyo+ghk1vpWmToeLT0kIGVMqZPk9viTmPj1g8F+97xD8BP7yyv1KdMe1sVnjJjOyqUmW1FM6PCiXeWI7xt+acCCthbZmqa3R8Qjo+PF3sI5SQ8hFjVcW/9V4q0TBY6A5CkZZCsBTieQXxKHpo4GS/ddM7JxMLV0yDX60z7aRyUeH/KKnBpuZUaw4sscEojh8bAQMIABh6QS3i8yKYhPj6pWl+jNSjrNVrCdU1F1fb7Yil/B+XJjJMmy5gZnda0n0bhiAWhUdCn5xV5LeXLdEUyY+lq470sKzM2U/R8fX0BAI8fP4a/v7+p/fHjx/Dx8SlxP0v8R8WsU8KL06s6bi45DTKdHrfS9DDwPBgJD44DUnQM5DoJOB6IjV2Mhsn3BPutm7wcj5u2RLCPvFILQqsAfNpEX+ZCupZyIWtNEvs1WkNmxP4dFVaezGwxk5mPQ5dDrTKTmXJmuryZsXSW8l7azK2FGjVqBF9fX5w4ccLUlpeXh9OnT6NDhw4ijswy+SpZ3MvhoDcYkMvx4AyAngd0PJCUa8Co/9uBnhd+EOzzda8ROP7Ca1VeG7BgId2VwR4Ib+NqdeG1FZSZiik7M5+h15/fC/b5utfb+PGFfpQZC2JVRS87OxuXL1/G5cuXYTAYcO/ePVy+fBl3794FwzAIDQ3FunXrcPjwYVy9ehVTpkyBk5MTBg8eLPbQLc4olRL5egOy9YDWAHAAGAA8gOCEi/jgwGpB/7uqZ/DhkHdxNUWLu1lafHIli64VsgKUmepTWmY6JFzCBweiBP3vqlph6WDKjKWxqo8LFy5cQL9+/Uw/R0REICIiAsOGDUNMTAxmzpyJ3NxczJ07F+np6QgKCsL+/fvpeqP/FD292k/JIiVfDx6AhDF+AvLITsPn2+ZAxj0JZ5bSBXOnrkWilgXPMHig4bHrVi7+StUjprMHfeq0YJSZqilPZtxz0vH59tlmMrMOCTopZcbCVGrBaVL7qrrgdNEzzziex7VULfINAHiAYRjk5Ouxb1Moel8/Jdj37UnrcbBlNzj8d4YaYLzGy00hwbCmyhq5VshSjv/XJHt4jVUl5u+ovJnZu3kK+lz7RbDviIkf40Cr7qVmxlIWka8tlvL3Th837IS5C1wdpBJI9RxkUhZaAzD32JZiBW9ttzE41LIbOB7I53gopcb9GYaBzgC6VojYrPJkJvy72GIFb13X0TjYqjtlxkJR0bMThS9w1XIG3MvmkKvjoOEYtHaToMP1c5h9eINgn9NPtUXEwDAYAOgMwsfjeR4yCUPXChGbVXZmzmPOofWCfc40bovlg2ZRZiyYVZ3IQirPV8mC43loOQOupOqRpuWhNTBwkwH59x7inTVhYPknKU129sDkCWvAyxRwZI3zFxxvnLTneR5ggMbOLN3hmdis0jKTd/8R3lkzU5CZFCd3TJ4QTZmxcPRNz06MUilxIVmH6xk6GACA5yFhGDRQGBC5ORxeGU/WgTQwDELHrMQDDz8wACQSCdxkBvg6slDKJQDP40U/Oaa0dKEJeWKzSspMfQcekVHFMzNlTCTu16lLmbFw9Nu3Ez5KKdZ0dMPwH1Og43nIJQzqO7F48/PVeCH+vKDvxtdCcb7Ni/CQGJdVUkgY1HWS4aW6DrTALbEbJWVmyK5ovHjznKDvJ30n42zblygzVoCKnhUraYX3kvgopejd0BFn1FqwDIPn/vwJU77ZJOhzoVVHHBg2A0wO0MhVajprraoX1xJiCaqamaALJzH1SIygz8VWwdg/fCZlxkpQ0bNSlbnHllqjR7bWgPh0PfzTH2Lux8L7tqW6eyMybC1YmRQDG7NwlrOmZY9ebaCo0D8WhFia6sjMO+YyM5MyY03oHbBSFb3HVuHAByg5REaEwS073bSdYyRYMmMNfmfcoE/Xo5W79Mm9wegmlsQGVCUzKqUBkStmwS0rzbSdYyT4cPoanJe4U2asCJ29aaUqeo+twoGf/L/VaHf7kmD7t2PDcdj/OYDn0cSVxd/pxtAWHA4q6R8LQqxFlTLz5Wo8+89FwfZvQ+bgUEPKjLWhjxxW6maGcaV3pezJNT+l3WMrScOBM/AI+PlbvH5km2CbrmdPnBsxBU2T9YJ/FApCSjexJNaq8BxeVTLzxtdbBdt0PXrg3MiplBkrRN/0rBQPHjczOGh0xhCVNXHuJGWQff0fRGydL2hP9amP3E2b8CiPLzGkBdcrFUY3sSSWruAQ4xm1Fv/mGCqVmawbt7F823uC9jTvesjdvJkyY6Wo6FkpRymLZu4sJAyDRs4SBPvIS50vkGrzEBs7C6552aY2LSvFZwvWg69Tp9SQjlIpoZAwpu10ZhqxBkUPMVYqM1tmwS03y9Smk0jx2YINlBkrRoc3rZijlEUjZwlWBnuU2fe1jR+i9d1rgrbNb7+L6wHtADy5ELfw4roFIS24Xsnab2JJ7Iu5Q4wVyUzfjUvR5u5VQdvmEe/gWjPKjDWjd8CKlfdwiSwuDl2OfiFo+yW4Nw6/NgbB/+1fVkgLbmJJiLXwVbK4ky0sfOXOzFdf4eWjuwVtpzr0wqHXQigzVo6KnpUq7+ESyY0bcJw1S9D2wK8hoiZHQMFKBPtTSIktKe2bWGkkN2/CMSxM0PbAtyFWha6gzNgAKnpWKthHbvZwSeGz1Tx0GsyaPgIuOTmm7Tq5AlsXxqD1U550uIXYtPIeYiyambAZxTOzbeEnlBkbQe+elSrrYlqd3oDQtfPhlxgv6JP6UQRmvP1CbQ2TEFGV9U2seGbeQ92Em4I+qcuWY/qIF2t6qKSW0NmbNqTw2Wqdju3B0LOHBdt/6DwAG55/Q6TREWJ5Cmem43dxGHr2kGD7D537Y0P7wSKNjtQEKno2pOBstSZ3rmLR7qWCbYkNVFg/aRnUeYYS9ibE/hRk5qmEa/igWGaaYsNEyoytoaJnQ3yVLBQ5mXh/9TQ46LWmdo3cER/NWQ+NwpEujiWkEEFmdPmmdmNmNiDHQUmZsTE0p2fFit4m5dX6cvSa9x7qP0oU9Fs7cSnuNGhKF8cSu2cuMz3nvY8GDxME/dZN+JAyY6Oo6Fkpc6u4198Zi8mnvxX0+63vMCS+9gaC6eJYYufMZ2YrJp8+Kuh3+tVhSOg3mDJjo+jdtFJFl1gKjL+EcTsiBH3+bRKI72cuQXgrFwousXtFM9Pi1mWM27Fc0OffJoH4vzDKjC2jOT0rVXiJJeesdLwXPR0yTmfanuPojA9mrcepDInpdieE2LOimXm/WGacKDN2gD7KWKmCJZakPI+568Phm/xAsP2dkOW4490QDgyDTC2HCT+nIcBNSndwJnarcGbCN8yF7+P7gu3vjqHM2AP6pmelClZxH3BgMzr8eUKwbVO3UYhr9QqupemRkafHjXQOtzN1+DfHgDNqLX2KJXbJlJmDnyL4j+OCbZu7jsSeZ3pQZuwAFT0rtepSFrrcOofxe9YI2s81bo2F/WdBo+eRqTXgUqoeWs4AB9b4VtMdnIm9UWv0iLqUiVWXsvDSP+cx/stowfbzjZ7BggGzKTN2gr6vW6nM+2qMWj4DUsOTOzGnOrlh7NjVSIMcjAGQMICeAzge8HZ4stI83cGZ2IvCZ2zWyUwplpk0pRvGjoumzNgRKnpWwNxhlblrZ8E347GgberoFbjjXg9SAyBljcFVsIADC6jzeLgojP3oDs7E1hVcj3fs3zykaQ1o6GjMjF+6WtBv6ugVuE2ZsStU9CxcwSfVooL+/k3w85qeE/DDM13AaXlIWcBRKgHAQ8kyeJxngDqXQ30nFjKWoQtuiU0r/O1OnctBo+fRe/8GPPeXMDNre47H960pM/aGip4FKLpKROEzxQquLSrNz02fR2T/GXCS8OBlDBQSQCFlkKXloeEARxbgASRmcRjYWIEpLekaJGLdypMZlmGgYBk8//evCD8aI9j/l6efw4r+MykzdohOZBFZwafSM2qt2TPFCl9bZHZ/F0+ET4qCQiGDt1KKr3t5op23AhLwAMMAPA+pRILWnnI0dZfCWc5SeIlVq0hm2mofY9O2dyDhn3xwVLt4InzyasqMnaKiJ7Kiq0QUPVPMV8mC481/0+MYCSJmrkGdpxogwF2K3g0dEeipwJqObvBQSKCUMnBXSNDCQwoHqYQm44lNKG9mWL0Oi9bNgld2qmlfA8NgxYxoeFBm7BZ9fBGZuW9yhYM2SqXEhWSduV2xc+hMXGv7AjieF8w5+CiNYT6j1goeu+hkfGmHiAixVOXNzMhd0Wh1/XdBv51vzsTVdi9SZuwYvVsiK1gloqSg+Sil+DTzN9THs4L9cl7uhjuTZqJRvrFv0fAVBL/gE3HRkJtbfPdCsg5rOrpRiIlFK29m6h36VLBfTpeuuD05jDJj5+idEllZQWMSElA3bBqwWHjmGbd1C+Z4uhd7vMKfRFWuLADjxHzRkJd2iCi8jWuNvmZCqqLMzCQmGjNTiKFePWNmvNyLPR5lxr5Q0ROZj1KKNR3dsDNeA3UuJwxafj6UISFgMopfssB7ehZrK/pJtOAfA3OfRMs6RESIpapoZniWhWbbNvBeXsUeizJjf6joWQAfpbTYJ0W1Ro+Uqe8g+MKFcj9ORT6JlnWIiBBLVnJm3kXwn38K2vMWLwYXHGz2cSgz9scmz96MjY1F69at4evriy5duuC3334reycLotbosT9qN4IPfFah/SrySbRg8d2CM0OLHiIi9sUWMrMv6gsEH9ghaM/u1QfaadPM7wTKjD2yuaK3f/9+zJs3D3PmzMHPP/+M9u3bY8iQIbh7967YQyu3oz/+hemfvFdqn6hLmZh7Og1RlzJN1yeZu7yhpE+iBYeIgn3kaOQsQbCPnCbk7ZQtZOab439jesx8QdtD7wZYNz0SYBjTotOUGcKkp6eXvtyHlenevTtatmyJjz/+2NT27LPPYsCAAfjggw9EHFk5aTTI6tQVDRJumJq0UjkUq4SHOV9poCg2BwGg3PMTli4+Ph4qlUrsYdQoS3mNlpyZcv2OcnOR1fHlYpmZs2wPtG3bIry1S4m5AMTLjPt24f380kPq1+jzic1S/t5tquhptVrUrVsXW7duxcCBA03t4eHhuHr1Ko4ePVpsn6J/eIRYkpr+h5AyQ2xNWZmxqcObKSkp4DgO3t7egnZvb2+o1eoS9iLEflFmiL2xqaJHCCGElMamip6npydYlsXjx8L7zD1+/Bg+Pj4ijYoQy0WZIfbGus5uKINcLkfbtm1x4sQJwfzEiRMn0L9/f7P7iD15HHUxA93en4quvx4RtEeERmLwwnHYGa/BGbUWP9zLF2zvUleOZu7SYqtGRF3KNLt+YLCP3KpWjbCUSe+aZAmv0dIzY+53VFJmVkxegTcWjTdlpnAGcvUcGDAWlRk6kUUcNlX0AGDq1KmYNGkSgoKC0KFDB2zbtg2PHj1CSEiI2EMza/pve+BTJLxfdHodx7oMwq+nM+AhZ8zeWqiZuxQrgz2KtZe1RBMhRVldZk4Xz8z/Or6Ob19+HadKyIyjlEUjZwllhthe0Xv99deRmpqKVatWISkpCS1atEBcXBwaNmwo9tCKYS9cgPfiBYK26w0C8PnkxXCUssg38LibY4CMLV70SloFotQlmggxw5oyI7l4Ed4fFM/MztDimSnvyimUGftiU5csWJX0dLi89BIk//5ratI4OGH6ykO4X+8pU5u3AkjTAt/8myfY/eZQX5sOpaUcCqlJ9vAaq0rwO0pPh0uXLpAkJpq2axycMCPyIO7Vb2JqK8iMpV+vSoc3xWFTJ7JYDZ6HcsoUQcEDgNWhywUFj+N5NHaVmS6iLcySwktIjeN5KKdOFRQ8AIievFxQ8ApnhlZOIebQX4EI5Bs2QFbkot/0kPE4++Jr4MzMK1BYib2Tb9wI2TffCNrSx4zDmc4lZ8aaTtwitYf+Na1l7JkzcFi8WNCmf/ZZMCs+whqOpXkFQopgz54tnpl27cBELqfMkAqjv45axCQnQzl2LBjuyQruvJsbNNu3AwoFfAD6dEpIIdK0NOP98fR6UxtlhlQFFb0aUPhOzL7K/z59OkjgOHEiJA8eCPpqYmLAN2ok0kgJsQwlZeapRYvMZ6ZxY3EGSqweFb1qVvROzHeyOVxI1uGzM1sgO35c0Dd/xgzoX31VpJESYhlKysyOs7HwOnNG0Dd/+nTKDKkSOnuzmpm7E3PgpV/hERUp6Kfv2BF5CxeKMURCLIq5zLS49BvqFM1McDDyFi0SY4jEhlDRq2ZF78RcJzUJ89fNgqTQjSoNXl7QbN0KyGRiDJEQi1IsM2lqzP94FiQGg6nN4OlJmSHVgopeNSt8J2YJp8e8tWHwyEgxbecZBrlbtoCvV0+sIRJiUYplZk0Y6qQnm7abMlPfti/eJrWDil41G6VSQiExXjM06su1aH31nGB7/jvvQN+1q0ijI8TyFM7MyD1r0frqWcH2/Llzoe/WTaTREVtDRa+aFazjN+qfX/DWgRjBNt3LLyP/nXdEGhkhlqkgMyNvn8Kw/cLMZD7/PPLffVekkRFbREWvBvgmP8DbK2YJ2gx+fsjdsgVgzS96S4g98015iBFmMnN76VLKDKlWVPSqm1YL5dixkKSnm5p4loVm2zbw3t7ijYsQS6XVQhkSAklamqmJZ1lotm6F3tNTxIERW0RFr5o5LFoE6e+/C9ryFi0C16mTSCMixLI5fPBB8cwsXAjuhRdEGhGxZVT0qpH00CEoNm0StOl69YJ2+nSRRkSIZZMePgxFTJG57169oJ0xQ6QREVtHRa+aSG7fhrJIcTP4+yN30yZAQr9mQoqS3L4N5bRpgjZDgwbIjYmhzJAaQ39Z1SE3F8rRo8FkZpqaeJkMmh07wHt4iDgwQixUXl7JmalTR8SBEVtHRa8aOM6bB/avvwRtjxd9CC4oSKQREWLZHMxlZuEScM89J9KIiL2goldFsi+/hPyzzwRtJzv2wdgWQ6DW6EvYixD7JduzB4odOwRtPwf3xtjANykzpMZR0asCybVrcJw9W9B2368R1oVGIJ83LqRLCHlCcv06HGcJr8e779cIa0NXUGZIraBbC1VWdjaUY8aA0TwJab5cgWXhG6FRuoAFoM7lSt6fEHuTk2OcxyuUGa1Mjo/mbIDGiTJDagd906sMnofj7Nlgb9wQNH8ybjHuNG4BAOB4Hj6OtJIEIQBKzcztpwIBUGZI7aCiVwmyzz6DPC5O0Pbjy4NwtOtgAMbwKiQMRqmUYgyPEIsj27kT8j17BG3HuwzEN92GAKDMkNpDhzcrSHLpEhyLLIDLtWiBRrFrEXyfgTqXg48ji1EqJXyU9OslRHL5MhyLLLTONW+OhlvXUWZIraO/sIrIyDDO4+Xnm5p4JydoduyAt5cbwr1EHBshligjwziPVzQzn31GmSGioMOb5cXzUE6bBvbOHUFz7tq1MDRrJtKgCLFgPA/l9OmUGWJRqOiVk3zTJsi+/lrQlj92LHRDhog0IkIsm3zzZsgOHxa05YeEUGaIqKjolQN7/jwcFi4UtHFt2iBv+XKRRkSIZWN//714Zlq3Rl5EhEgjIsSIil4ZmNRUKENCwOifrBTBu7oi57PPAAcHEUdGiGViUlONc986namNd3WFhjJDLAAVvdIYDHCcNAmSe/cEzZqNG8E3bizOmAixZAYDHCdPLp6ZDRtgeOopkQZFyBNU9EqhWLsWsu+/F7TlT50Kfb9+Io2IEMumWLcOsv/7P0Fb/pQp0PfvL9KICBGiolcC9tQpKJYtE7Tp27dH3uLF4gyIEAvHnjoFxdKlgjb9889TZohFoaJnBpOUBOW4cWAMBlOboU4daLZtA2QyEUdGiGVi1Goox48XZsbDw5gZuVzEkREiREWvKI6Dcvx4SJKSBM25n34KvkEDkQZFiAUryMyjR4Lm3E8/Be/vL9KgCDGPil4RirVrIf3lF0FbXng49K+8ItKICLFsinXrIP35Z0Fb3pw50PfoIdKICCkZFb0itMOGQd+xo+lnfefOyJ8/X8QREWLZtG+9JczMiy9SZojFoqJXBF+vHnIOH0b+zJkw+PlBExsLsHS7E0JKwterh5yvv0ZeWBgMvr7GzEhpWV9imajomSOTIW/JEmSdOQPe11fs0RBi+aRS5C9ejKyzZ8H7+Yk9GkJKREWvNO7uYo+AEOtCmSEWzmqK3o4dO/Daa6+hYcOGcHd3R2JiYrE+6enpmDhxIho2bIiGDRti4sSJSE9Pr/3BEmIBKDOEFGc1RU+j0aBbt26YN29eiX3Gjx+Py5cvY+/evdi7dy8uX76MSZMm1eIoCbEclBlCirOa2eYpU6YAAC5cuGB2+40bN/DDDz/g2LFjaN++PQBgzZo16NOnD+Lj46FSqcr1PGqNHjvjNUjScPBV0t2cifWizBBSnNV80yvLuXPn4OzsjA4dOpjagoOD4eTkhLNnz5brMdQaPWadzsAZtRb/5hhwRq3FrNMZUGv0Ze9MiJWhzBB7ZDMfx9RqNTw9PcEwjKmNYRh4eXlBrVaXuF98fLzp/7f+K0VaJgvJk4eAhgeiT2diXEOxQ6wU/FR43LaKXiPK/W2rMqwlM7b7d0CZrgllZUbUords2TJERUWV2ufrr79G586da2wMhX9BXHIaXHhDsT4GJwlUKo8aG0O5nLov+LEm/zG0BBU5vGatKvMa7S0zNv13QJkWhahFLzQ0FG+++WapfRqUc71LHx8fpKSkgOd50ydXnueRnJwMHx+fcj2Gr5LFnWwObKFPvhzPw8eRLk4nloEyQ0jViFr0PD094enpWS2P1b59e2RnZ+PcuXOmOYpz584hJydHMGdRmlEqJS4k65Bv4MEyDDieh0LCYJRKWfbOhNQCygwhVWM1c3pJSUlISkrCrVu3ABjPPMvIyIC/vz88PDzQrFkzvPLKK5g1axbWrl0LAJg1axZ69epV7q/UPkop1nR0w854DdS5HHwc6Uw0Yr0oM4QUZzV/mdu2bUNkZKTp54JDPBs3bsTbb78NAIiNjcU777yDN954AwDQp08frFy5skLP46OUIryNazWNmhDxUGYIKY5JT0/nxR4EKZv7duGkd3pIfZFGUjssZdK7JtnDa6wqW/4dUabFYTPX6RFCCCFloaJHCCHEblDRI4QQYjdoTo8QQojdoG96hBBC7AYVPUIIIXaDih4hhBC7QUWPEEKI3aCiRwghxG5Q0bMisbGxaN26NXx9fdGlSxf89ttvYg+p2kRERMDd3V3wX0BAgNjDqrJff/0Vb731Flq0aAF3d3fs3r1bsJ3neURERKB58+bw8/ND3759ce3aNZFGa7n69u1b7O9j7NixYg+ryijTtY+KnpXYv38/5s2bhzlz5uDnn39G+/btMWTIENy9e1fsoVUblUqFGzdumP6zhX8AcnJyEBgYiBUrVsDR0bHY9nXr1mHjxo2IjIzE8ePH4e3tjUGDBiErK0uE0Vq2t99+W/D3sWbNGrGHVCWUaXFQ0bMSGzduxPDhwzF69Gg0a9YMq1atgq+vL7Zt2yb20KqNVCqFr6+v6T8vLy+xh1RlPXv2xKJFizBgwABIJMK48TyPmJgYhIWFYcCAAQgMDERMTAyys7Oxd+9ekUZsuZRKpeDvw83NTewhVQllWhxU9KyAVqvFxYsX0a1bN0F7t27dcPbsWZFGVf0SEhLQvHlztG7dGmPHjkVCQoLYQ6pRiYmJSEpKEryvjo6O6NSpk029r9Vl3759aNKkCYKDg7FgwQKr/jZMmRaP1dxayJ6lpKSA4zh4e3sL2r29vaFWq0UaVfV67rnn8Mknn0ClUiE5ORmrVq1Cz549cebMGdSpU0fs4dWIpKQkADD7vj58+FCMIVmsIUOGwN/fH35+frh+/TqWLFmCK1eu4MCBA2IPrVIo0+JlmooesQg9evQQ/Pzcc8+hbdu2+OKLLzBt2jSRRkVq0rJlyxAVFVVqn6+//hqdO3fGmDFjTG0tW7ZE48aN0b17d1y8eBFt27at2YGSSrHUTFPRswKenp5gWRaPHz8WtD9+/Bg+Pj4ijapmOTs7o3nz5rh9+7bYQ6kxvr6+AIzvo7+/v6ndlt/XwkJDQ003ti1JgwYNzLa3a9cOLMvi9u3bVln0KNPioTk9KyCXy9G2bVucOHFC0H7ixAl06NBBpFHVrLy8PMTHx5sKgy1q1KgRfH19Be9rXl4eTp8+bbPva2Genp4ICAgo9T+lUml23ytXroDjOKv9+6BMi4e+6VmJqVOnYtKkSQgKCkKHDh2wbds2PHr0CCEhIWIPrVosWLAAvXv3RoMGDUzH/zUaDYYNGyb20KokOzvb9MnWYDDg3r17uHz5Mjw8PODv74/Q0FBER0dDpVKhadOmiIqKgpOTEwYPHizyyC3HnTt3EBcXh549e6JOnTq4ceMGFixYgNatWyM4OFjs4VUaZVocVPSsxOuvv47U1FSsWrUKSUlJaNGiBeLi4tCwYUOxh1YtHjx4gPHjxyMlJQVeXl547rnn8P3331v967tw4QL69etn+jkiIgIREREYNmwYYmJiMHPmTOTm5mLu3LlIT09HUFAQ9u/fDxcXFxFHbVlkMhlOnjyJTZs2IScnB/Xr10fPnj0xb948sCwr9vAqjTItDrqfHiGEELtBc3qEEELsBhU9QgghdoOKHiGEELtBRY8QQojdoKJHCCHEblDRI4QQYjeo6BFCCLEbVPQIITZl9+7d+OWXX8QeBrFQdHE6IcQmfPHFF1AqlcjJyYG/vz9u3LiBp556Cq+88orYQyMWhL7pEUJswptvvon79+8jJiYGS5cuhVKpRPfu3QV9Dh06BH9/f6Snp4szSCvUvHnzCt0KaMuWLWjVqhXy8/NrcFSVR0WPVNqDBw/g5+eH8ePHC9rPnTsHPz8/ug8eqXUMw5j+n2VZwc8cx2H58uUYO3Ys3N3dRRideStXrsTJkyfFHoZZ6enpePToEVq2bFnufUaMGIH8/Hxs3769BkdWeVT0SKXVq1cPo0ePxoEDB/DPP/8AABISEjB8+HC0b98ea9asEXmExJ7ExcWhbt26CA0NxcKFC5GZmYkffvjBtP3YsWO4ceOGRd3F4M6dO1i+fDnu3bsn9lDMunbtGgAgMDCw3Ps4Ojpi2LBh2LBhAwwGQ00NrdKo6JEqmT17NuRyOaKjo5GRkYGhQ4eiTp062LlzJ2QymdjDI3Zk+PDhGDRoEADjN76JEycK5vN27dqFdu3aoXHjxiKNsLiLFy8CQLXfCLe6Di1ev34dACr0TQ8w3kHi3r17+Pnnn6tlHNWJih6pEl9fX4wdOxZxcXF48803kZKSgri4OIs6fESsW0UPo7/99tvo3LmzoC0vLw8//vgjunbtWuLz/PXXXwgJCYFKpYKvry+CgoLwwQcfCPqcPn0aQ4cORaNGjdCwYUO8/vrrpsJQYODAgXjllVfw999/Y+jQofD390dAQACioqIE/Xr27Gn61vnCCy/A3d0d3t7e0Ov15X4uABg8eDC6dOmC33//Hf3790fdunUxderUCj0GABw+fBhdunSBn58fOnXqhF9++QXXrl2Dj48PvLy8KvQ+tG3bFh4eHjhy5EiJv2+x0NmbpMoePXqEwMBASKVSHDlyBO3btxd7SMTGvPvuu4iNjcW5c+fw9NNPIyEhAa+88goCAwOxb9++Mo8qnDlzBr1798bOnTvRv3//YtuPHj2KkJAQNGjQAKNGjYKLiwsuX76Ma9eu4bvvvgMA7NmzB6GhoejcuTP69OkDnU6H2NhYaDQanD9/3vRBLyAgAF5eXkhJScHw4cPRsGFDfPnllzh79iy++eYbvPDCCwCMRSYyMhJZWVl4//33AQBKpRL9+vUr93MBQIsWLaBUKpGSkoKRI0fi6aefRr169ZCWllbux4iNjUV4eDj69OmDXr164ebNm9i9ezfq168PHx8fHDx4sMLvw8CBA5GcnIxTp05V9O2uUXQTWVJlq1evhsFgAMdx8PDwEHs4xAbNnj0bO3fuRHR0NJYvX17hw+g3b94EALOHNm/fvo0JEyagW7du2LZtGxwdHU3bCg4T3rhxAzNmzMDs2bOxYMEC0/Z+/fqhXbt2OHjwIMaMGYPk5GSo1WpotVr89NNPpufr0aMHWrVqhQsXLpiKXv/+/bFo0SK0b98eQ4cONT1meZ8LANLS0vDw4UM4Ojri+++/R6tWrUyPMXLkyHI9xtWrVzF//nxMmzYNy5YtM/WVSqX4+OOP8fLLL1fqfWjcuDF+++23Ut8XMdDhTVIln3zyCbZs2YIVK1bAyckJkZGRYg+J2KCqHkZPS0sDALP9V6xYAalUipiYGEHBAwCFQgHAeIalt7c3Jk6ciJSUFNN/Li4u8PDwQEJCAgDgypUrAIDw8HBBgS0oCEql0tSWnZ2NxMREU6EqUN7nAoC///4bADBz5kzB41TkMVavXg2lUol3331XMI6C4lz4JJaKvA/u7u7QarXIzMwstk1M9E2PVNrRo0exYMEChIeHY/LkyUhOTkZ0dDTmzp2LZs2aiT08YmOmT5+OmJgYXLhwAUeOHKnUCSk8L5zNyc/PxzfffINhw4aVWEDz8/Nx7Ngx5OTkICAgwGwfFxcXAMZvTYDxG1VhBd80mzZtamq7evUqeJ4XFKuKPFfh53vjjTcq9RharRbHjh3D8OHDBY8LwDS3WLQol/d9KPhdF75sxBJQ0SOVcvHiRUyYMAEDBgwwzUdMmTIFmzZtwooVKyz2Gh1ivapyGL1OnToAjNedNWrUyNSekJCAnJwctGnTpsR9ExMTkZOTg3nz5iE4ONhsn+bNmwMwftPz9PQUPAfw5BtZ4QJS8K2wcFtFnqvgMby8vKBSqSr1GHfu3EFOTo7Zs0cvXLgAiURS7ANsed+H9PR0yOXyYsVUbFT0SIXdv38fw4YNQ2BgIGJiYkyf5OrUqYOxY8di/fr1mDt3boWu7SGkNIUPoy9fvhyRkZGIjY0t9/4F/3AnJiYKClxubi6A0r+NZGVlAQBUKpVgfsucq1evFvtmBBiLU926dU3Ft6DN09MTfn5+lXqugsco+nwVeYyC11+UVqvFrl270KRJE8Eh34q8D4mJiSV+0xQTzemRCqtfvz6uXbuG77//Hg4ODoJtH374IdLS0qjgkWpT9DD6xIkTsX//fty4caPcj9G2bVs4ODiYrosr0KRJEzAMY3ZFFJ1OBwBo2LAhGIbBoUOHivUxGAxITU01/f/169fNFr2///672LVud+/eRf369QVt5X0uwHj40NzzVeQxGjRoAADFrqdbt26d6azsAhV9Hy5duoQOHTqY3SYm+qZHCLFY1XUYXS6Xo3v37jhx4gQWLVpkand1dcWIESPw+eefQ6fT4aWXXkJOTg5+++03PPfcc5gzZw68vb0xZMgQxMXF4Y033kCvXr3A8zwSExNx5MgRrFy5Er1790ZCQgI0Gk2xIsRxHK5fv45JkyYJ2hs1aoTjx49jzZo1qFevHpo2bYqgoKByPRfw5NBk0ecr73gBwMvLCz169EBcXBwcHBzw7LPP4uTJk6YiWFD0Kvo+XLx4EWlpaejbt2+53p/aREWPEGKRqvsw+ogRIzBs2DAkJCQITr5YuXIl6tWrh3379uHYsWNwdXVFUFAQ+vTpY+qzYcMGtGzZEl9++SUWL14MhUKBRo0a4c0330THjh0BmJ+jA4Bbt24hLy+v2De9OXPm4Pbt24iOjkZWVhYWL16MoKCgcj1X4eczt1pKeR8DMB6yDAsLw4EDB3Do0CF07twZq1evxujRoxEYGFip9+HAgQNo0KBBuQ7R1ja6OJ0QYhcMBgM6duyI3r17Y8mSJWIPx2bl5eWhdevWCAsLw5QpU8QeTjE0p0cIsQsSiQTvvfcetm3bRrcWqkGff/455HI5xo0bJ/ZQzKJveoQQQuwGfdMjhBBiN6joEUIIsRtU9AghhNgNKnqEEELsBhU9QgghdoOKHiGEELtBRY8QQojdoKJHCCHEblDRI4QQYjf+H3E1dUmqcW5tAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x     = np.random.normal(3,3,100) # generate 1000 points\n",
    "xbar  = np.mean(x)\n",
    "xstar = x-xbar\n",
    "\n",
    "epsilons = np.random.normal(0,2,100)\n",
    "\n",
    "y= 2+3*x + epsilons\n",
    "\n",
    "fig,axs =plt.subplots(1,2) \n",
    "ax=axs[0]\n",
    "ax.scatter(x,y,alpha=0.7)\n",
    "ax.set(xlabel=r\"$x$\",ylabel=\"y\")\n",
    "\n",
    "minx,maxx = min(x),max(x)\n",
    "b1,b0 = np.polyfit(x,y,1)\n",
    "ax.plot([minx,maxx],[b0+b1*minx,b0+b1*maxx],color=\"red\" )\n",
    "\n",
    "ax.axvline(0)\n",
    "ax.axhline(0)\n",
    "\n",
    "ax=axs[1]\n",
    "ax.scatter(xstar,y,alpha=0.7)\n",
    "ax.set(xlabel=r\"$x^{*} (centered x)$\",ylabel=\"y\")\n",
    "\n",
    "minx,maxx = min(xstar),max(xstar)\n",
    "b1,b0 = np.polyfit(xstar,y,1)\n",
    "ax.plot([minx,maxx],[b0+b1*minx,b0+b1*maxx],color=\"red\" )\n",
    "\n",
    "ax.axvline(0)\n",
    "ax.axhline(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Are $\\beta_{0}$ and $\\beta_{0}^{*}$ equivalent?\n",
    "\n",
    "The intercept for the SLR and SLR with centered $x$ are **different**.\n",
    "\n",
    "From above, we saw that the SLR with centered x simplifed to\n",
    "###### Algebraic\n",
    "\\begin{align}\n",
    "    Y_{i} &=  \\beta^{*}_{0}-\\beta^{*}_{1}\\bar{x} + \\beta^{*}_{1}x_{i} + \\epsilon_{i}\\\\\n",
    "    \\epsilon_{i} &\\sim \\mathcal{N}\\left(0,\\sigma^{2}\\right)\\\\\n",
    "\\end{align}\n",
    "\n",
    "The expression $\\beta^{*}_{0}-\\beta^{*}_{1}\\bar{x}$ involve only constants.\n",
    "This expression is $\\beta_{0}$. That is, $\\beta_{0} = \\beta^{*}_{0}-\\beta^{*}_{1}\\bar{x}$\n",
    "\n",
    "###### Visually\n",
    "We see from above that the estimated regression line (red line) intersects the y-axis at different points depending on whether we center or do not center our data.\n",
    "\n",
    "Though centering makes computing MLEs easier this process also changes our $\\beta_{0}$ or y-intercept."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix algebra to compute the MLEs for centered SLR\n",
    "\n",
    "Let us look again at the MLE estimates for sentered SLR\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{\\beta_{0}^{*}} &= \\bar{y} \\\\\n",
    "    \\hat{\\beta_{1}^{*}} &= \\frac{ \\sum_{i=1}^{N}x^{*}_{i}y_{i}}{\\sum_{i=1}^{N}(x^{*}_{i})^{2}}\n",
    "\\end{align}\n",
    "\n",
    "The MLE for $\\hat{\\beta_{1}^{*}}$ appears to be two inner products divided by one another. \n",
    "That is, we can compute $\\hat{\\beta_{1}^{*}}$ as \n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{\\beta_{1}^{*}} = \\frac{x^{*} \\cdot y}{x^{*} \\cdot x^{*} }\n",
    "\\end{align}\n",
    "\n",
    "where $x \\cdot y = \\sum_{i=1}^{N}x_{i}y_{i}$ and \n",
    "\n",
    "\\begin{align}\n",
    "        x = \\begin{bmatrix}\n",
    "                        x_{1}\\\\\n",
    "                        x_{2}\\\\\n",
    "                        x_{3}\\\\\n",
    "                        \\vdots\\\\\n",
    "                        x_{N}\n",
    "                      \\end{bmatrix} \\; \\; \n",
    "        y = \\begin{bmatrix}\n",
    "                        y_{1}\\\\\n",
    "                        y_{2}\\\\\n",
    "                        y_{3}\\\\\n",
    "                        \\vdots\\\\\n",
    "                        y_{N}\n",
    "                      \\end{bmatrix}\n",
    "\\end{align}              \n",
    "\n",
    "We can also recognize the MLE for $\\hat{\\beta^{*}_{0}}$ as an inner product\n",
    "\\begin{align}\n",
    "    \\hat{\\beta_{0}^{*}} &= \\bar{y} \\\\\n",
    "                        &= \\frac{1}{N} \\sum_{i=1}^{N} y_{i}\\\\\n",
    "                        &= \\frac{1}{N} \\sum_{i=1}^{N} 1 \\times y_{i}\\\\\n",
    "                        &= \\frac{1}{N} \\mathbb{1} \\cdot y\n",
    "\\end{align}\n",
    "where \n",
    "\n",
    "\\begin{align}\n",
    "        \\mathbb{1} = \\begin{bmatrix}\n",
    "                        1\\\\\n",
    "                        1\\\\\n",
    "                        1\\\\\n",
    "                        \\vdots\\\\\n",
    "                        1\n",
    "                      \\end{bmatrix}\n",
    "\\end{align}   \n",
    "\n",
    "Using matrix algebra, our MLEs become\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{\\beta_{0}^{*}} &= (\\mathbb{1} \\cdot \\mathbb{1})^{-1} (\\mathbb{1} \\cdot y) \\\\ \n",
    "    \\hat{\\beta_{1}^{*}} &= (x^{*} \\cdot x^{*})^{-1} (x^{*} \\cdot y )\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A data setup such that we arrive at the above MLEs\n",
    "\n",
    "The above MLE estimates pop out of the following data setup. \n",
    "Let $X$ be a matrix with $N$ rows and $2$ columns. \n",
    "The first column will be all ones and the second column will be filled, from top to bottom, with $x_{1},x_{2},\\cdots,x_{N}$. \n",
    "\n",
    "\\begin{equation}\n",
    "    X = \\begin{bmatrix}\n",
    "            1 & x_{1} \\\\ \n",
    "            1 & x_{2} \\\\ \n",
    "            1 & x_{3} \\\\ \n",
    "            \\vdots & \\vdots \\\\ \n",
    "            1 & x_{N} \\\\ \n",
    "        \\end{bmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "The above MLE formulas may motivate us to try \n",
    "\n",
    "\\begin{align}\n",
    "    (X'X)^{-1} X'y\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "    X'y = \\begin{bmatrix}\n",
    "            1     & 1 & 1  & \\cdots  &1 \\\\\n",
    "            x_{1} & x_{2} & x_{3} & \\cdots &x_{N}\\\\\n",
    "        \\end{bmatrix} \\times\n",
    "        \\begin{bmatrix}\n",
    "                        y_{1}\\\\\n",
    "                        y_{2}\\\\\n",
    "                        y_{3}\\\\\n",
    "                        \\vdots\\\\\n",
    "                        y_{N}\n",
    "                      \\end{bmatrix} = \n",
    "                      \\begin{bmatrix}\n",
    "                          \\mathbb{1}' y \\\\\n",
    "                          x'y\n",
    "                      \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "The first multiplication looks like a step in the right direction. \n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "    X'X = \\begin{bmatrix}\n",
    "            1     & 1 & 1  & \\cdots  &1 \\\\\n",
    "            x_{1} & x_{2} & x_{3} & \\cdots &x_{N}\\\\\n",
    "        \\end{bmatrix} \\times\n",
    "        \\begin{bmatrix}\n",
    "            1 & x_{1} \\\\ \n",
    "            1 & x_{2} \\\\ \n",
    "            1 & x_{3} \\\\ \n",
    "            \\vdots & \\vdots \\\\ \n",
    "            1 & x_{N} \\\\ \n",
    "        \\end{bmatrix} = \n",
    "        \\begin{bmatrix}\n",
    "            \\mathbb{1}'\\mathbb{1} & \\mathbb{1}'x\\\\\n",
    "             \\mathbb{1}'x         & x'x\\\\\n",
    "        \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "If we centered our x covariate then $\\bar{x}$ is zero and so \n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbb{1}'x = \\sum_{i=1}^{N} x_{i} =N\\bar{x} = 0\n",
    "\\end{align}\n",
    "\n",
    "We can use the above to simplify $X'X$\n",
    "\n",
    "\\begin{align}\n",
    "    X'X = \n",
    "        \\begin{bmatrix}\n",
    "            \\mathbb{1}'\\mathbb{1} & 0\\\\\n",
    "             0         & x'x\\\\\n",
    "        \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "The matrix $X'X$ is now a diagonal matrix---a matrix with non-zero elements for the entries $(i,i)$.\n",
    "The inverse of a diagonal matrix is easy to compute\n",
    "\n",
    "\\begin{align}\n",
    "    (X'X)^{-1} = \n",
    "        \\begin{bmatrix}\n",
    "            1 \\Big/ \\mathbb{1}'\\mathbb{1} & 0\\\\\n",
    "             0         & 1 \\Big/ x'x\\\\\n",
    "        \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "and we can now see our final product is \n",
    "\n",
    "\\begin{align}\n",
    "    (X'X)^{-1}X'y =\n",
    "        \\begin{bmatrix}\n",
    "            1 \\Big/ \\mathbb{1}'\\mathbb{1} & 0\\\\\n",
    "             0         & 1 \\Big/ x'x\\\\\n",
    "        \\end{bmatrix}\n",
    "           %\n",
    "        \\begin{bmatrix}\n",
    "            \\mathbb{1}' y \\\\\n",
    "            x'y\n",
    "         \\end{bmatrix} = \n",
    "         \\begin{bmatrix}\n",
    "            \\mathbb{1}' y \\Big/ \\mathbb{1}'\\mathbb{1} \\\\ \n",
    "            x'y\\Big/ x'x\n",
    "         \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "We can see that the above matrix multiplication returns a vector with two MLEs: $\\hat{\\beta^{*}_{0}}$ and $\\hat{\\beta^{*}_{1}}$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Variance for the MLEs follows a similar pattern. \n",
    "\n",
    "We saw in simple linear regression that the variance for $\\hat{\\beta_{1}}$ was \n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathbb{V}(\\hat{\\beta_{1}})= \\sigma^{2}\\frac{1}{\\sum_{i=1}^{N}(x_{i}-\\bar{x})^{2}}\n",
    "\\end{equation}\n",
    "\n",
    "and the variance for $\\hat{\\beta_{0}}$\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathbb{V}(\\hat{\\beta_{0}})= \\sigma^{2}\\left( \\frac{1}{N} + \\frac{1}{\\sum_{i=1}^{N}(x_{i}-\\bar{x})^{2}} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "There is a natural generaliztion for the variance of $\\hat{\\beta_{0}}$ and $\\hat{\\beta_{1}}$ like the above generaliztion for the MLEs. \n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathbb{V}(\\hat{\\beta}) = \\sigma^{2}(X'X)^{-1}\n",
    "\\end{equation}\n",
    "\n",
    "The above variance will return a **covariance matrix** where the variance of $\\hat{beta_{0}}$ is the first diagonal entry and the variance of $\\hat{\\beta_{1}}$ is the second diagonal entry. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Linear Regression (finally) and an example\n",
    "\n",
    "Suppose we collect data points\n",
    "\\begin{equation}\n",
    "    \\mathcal{D} = [ (x^{1}_{1},x^{2}_{1},x^{3}_{1},y_{1}),(x^{1}_{2},x^{2}_{2},x^{3}_{2},y_{2}),\\cdots,(x^{1}_{N},x^{2}_{N},x^{3}_{N},y_{N})  ]\n",
    "\\end{equation}\n",
    "\n",
    "and propose that the conditional distribution of $Y_{i}$ is\n",
    "\\begin{align}\n",
    "    Y_{i} | x_{i}^{1},x_{i}^{2},\\cdots,x_{i}^{3} \\sim \\mathcal{N}\\left( \\beta_{0} + \\beta_{1}x_{i}^{1}+ \\beta_{2}x_{i}^{2}+ \\beta_{3}x_{i}^{3}, \\sigma^{2} \\right)\n",
    "\\end{align}\n",
    "\n",
    "The above is **multivariate linear regression** model. \n",
    "We can compute the MLE for $\\hat{\\beta}$ as \n",
    "\n",
    "\\begin{equation}\n",
    "    \\hat{\\beta} = (X'X)^{-1}X'y\n",
    "\\end{equation}\n",
    "\n",
    "and we can compute the variance for our MLEs as \n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathbb{V}(\\hat{\\beta}) = \\sigma^{2}(X'X)^{-1}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_diabetes \n",
    "import numpy.linalg as la\n",
    "\n",
    "\n",
    "diabData = load_diabetes()\n",
    "X = diabData[\"data\"]\n",
    "y = diabData[\"target\"]\n",
    "\n",
    "# add intercept (column of ones)\n",
    "X = np.column_stack([np.ones((len(X),1)),X]  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "betaHat = la.solve( X.T.dot(X), X.T.dot(y))\n",
    "betaHat\n",
    "\n",
    "epsilons = y - X.dot(betaHat)\n",
    "sigma2Hat = epsilons.T.dot(epsilons) / len(epsilons)\n",
    "\n",
    "varBetaHat = sigma2Hat * la.inv(X.T.dot(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "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.518</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.507</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   46.27</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 26 Oct 2021</td> <th>  Prob (F-statistic):</th> <td>3.83e-62</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>00:27:29</td>     <th>  Log-Likelihood:    </th> <td> -2386.0</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   442</td>      <th>  AIC:               </th> <td>   4794.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   431</td>      <th>  BIC:               </th> <td>   4839.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    10</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>  152.1335</td> <td>    2.576</td> <td>   59.061</td> <td> 0.000</td> <td>  147.071</td> <td>  157.196</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>  -10.0122</td> <td>   59.749</td> <td>   -0.168</td> <td> 0.867</td> <td> -127.448</td> <td>  107.424</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>    <td> -239.8191</td> <td>   61.222</td> <td>   -3.917</td> <td> 0.000</td> <td> -360.151</td> <td> -119.488</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x3</th>    <td>  519.8398</td> <td>   66.534</td> <td>    7.813</td> <td> 0.000</td> <td>  389.069</td> <td>  650.610</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x4</th>    <td>  324.3904</td> <td>   65.422</td> <td>    4.958</td> <td> 0.000</td> <td>  195.805</td> <td>  452.976</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x5</th>    <td> -792.1842</td> <td>  416.684</td> <td>   -1.901</td> <td> 0.058</td> <td>-1611.169</td> <td>   26.801</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x6</th>    <td>  476.7458</td> <td>  339.035</td> <td>    1.406</td> <td> 0.160</td> <td> -189.621</td> <td> 1143.113</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x7</th>    <td>  101.0446</td> <td>  212.533</td> <td>    0.475</td> <td> 0.635</td> <td> -316.685</td> <td>  518.774</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x8</th>    <td>  177.0642</td> <td>  161.476</td> <td>    1.097</td> <td> 0.273</td> <td> -140.313</td> <td>  494.442</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x9</th>    <td>  751.2793</td> <td>  171.902</td> <td>    4.370</td> <td> 0.000</td> <td>  413.409</td> <td> 1089.150</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x10</th>   <td>   67.6254</td> <td>   65.984</td> <td>    1.025</td> <td> 0.306</td> <td>  -62.065</td> <td>  197.316</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.506</td> <th>  Durbin-Watson:     </th> <td>   2.029</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.471</td> <th>  Jarque-Bera (JB):  </th> <td>   1.404</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.017</td> <th>  Prob(JB):          </th> <td>   0.496</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.726</td> <th>  Cond. No.          </th> <td>    227.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<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.518\n",
       "Model:                            OLS   Adj. R-squared:                  0.507\n",
       "Method:                 Least Squares   F-statistic:                     46.27\n",
       "Date:                Tue, 26 Oct 2021   Prob (F-statistic):           3.83e-62\n",
       "Time:                        00:27:29   Log-Likelihood:                -2386.0\n",
       "No. Observations:                 442   AIC:                             4794.\n",
       "Df Residuals:                     431   BIC:                             4839.\n",
       "Df Model:                          10                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const        152.1335      2.576     59.061      0.000     147.071     157.196\n",
       "x1           -10.0122     59.749     -0.168      0.867    -127.448     107.424\n",
       "x2          -239.8191     61.222     -3.917      0.000    -360.151    -119.488\n",
       "x3           519.8398     66.534      7.813      0.000     389.069     650.610\n",
       "x4           324.3904     65.422      4.958      0.000     195.805     452.976\n",
       "x5          -792.1842    416.684     -1.901      0.058   -1611.169      26.801\n",
       "x6           476.7458    339.035      1.406      0.160    -189.621    1143.113\n",
       "x7           101.0446    212.533      0.475      0.635    -316.685     518.774\n",
       "x8           177.0642    161.476      1.097      0.273    -140.313     494.442\n",
       "x9           751.2793    171.902      4.370      0.000     413.409    1089.150\n",
       "x10           67.6254     65.984      1.025      0.306     -62.065     197.316\n",
       "==============================================================================\n",
       "Omnibus:                        1.506   Durbin-Watson:                   2.029\n",
       "Prob(Omnibus):                  0.471   Jarque-Bera (JB):                1.404\n",
       "Skew:                           0.017   Prob(JB):                        0.496\n",
       "Kurtosis:                       2.726   Cond. No.                         227.\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "model = sm.OLS(y,X)\n",
    "res = model.fit().summary()\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-10.012197817471035\n",
      "59.00101633853383\n"
     ]
    }
   ],
   "source": [
    "print(betaHat[1])\n",
    "print(np.sqrt(varBetaHat[1,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}