chapter11.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup and imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 11.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 11.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = pd.DataFrame({\n",
    "    'A': [1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],\n",
    "    'Y': [\n",
    "        200, 150, 220, 110, 50, 180, 90, 170,\n",
    "        170, 30, 70, 110, 80, 50, 10, 20\n",
    "    ]\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAETRJREFUeJzt3X+s3XV9x/Hn21JdN9kK9GLgUq0uhck0s+TEsLA5DIsF/rDVgKGLoxqyLg6X/TBNYCbTbFkk69TEzOhqJBQzUTaxNJOtuurCtgy209XxQ9PZKUJvCb2KZWbcOajv/XG+Fy7lc889bc/3fM+53+cjuTnnfM733u/703s5L76f7/fz/URmIknSiV7SdAGSpPFkQEiSigwISVKRASFJKjIgJElFBoQkqciAkCQVGRCSpCIDQpJUdEbTBZyONWvW5Lp165ouQ5Imyv79+7+XmVNLbTfRAbFu3Tq63W7TZUjSRImI7w6ynUNMkqQiA0KSVGRASJKKDAhJUpEBIUkqmuirmCSpDXYfmGHH3oMcOTbH+atXsX3jRWzeMF37fg0ISRpjuw/McPNdDzL3zHEAZo7NcfNdDwLUHhIOMUnSGNux9+Bz4TBv7pnj7Nh7sPZ9GxCSNMaOHJs7qfZhMiAkaYydv3rVSbUPkwEhSWNs+8aLWLVyxQvaVq1cwfaNF9W+b09SS9IYmz8R7VVMkqQX2bxheiSBcCKHmCRJRQaEJKnIgJAkFRkQkqQiA0KSVGRASJKKDAhJUpEBIUkqMiAkSUUGhCSpqLaAiIi1EfG1iPhmRDwcEb9TtZ8dEV+JiG9Vj2dV7RERH4uIQxHxQERcUldtkqSl1XkE8Szwvsx8LXApcGNEXAzcBOzLzPXAvuo1wFXA+uprG/CJGmuTpImx+8AMl93yVV5905e47JavsvvAzEj2W1tAZObjmfnv1fMfAt8EpoFNwK5qs13A5ur5JuD27LkPWB0R59VVnyRNgvklR2eOzZE8v+ToKEJiJOcgImIdsAG4H3hFZj4OvRABzq02mwYeW/Bth6s2SWqtZb3kaES8HPgC8LuZ+d/9Ni20ZeHnbYuIbkR0Z2dnh1WmJI2lZbvkaESspBcOf5mZd1XNT8wPHVWPR6v2w8DaBd9+AXDkxJ+ZmTszs5OZnampqfqKl6QxsCyXHI2IAD4NfDMzP7LgrT3A1ur5VuDuBe3XV1czXQo8NT8UJUlttVyXHL0M+HXgwYj4etX2B8AtwJ0RcQPwKHBt9d49wNXAIeBp4N011iZJE6HJJUcj80XD/BOj0+lkt9ttugxJmigRsT8zO0tt50xqSVKRASFJKjIgJElFBoQkqciAkCQVGRCSpCIDQpJUZEBIkooMCElSkQEhSSoyICRJRXXerE+SNAS7D8w0crM+A0KSxtj8kqPzq8rNLzkK1B4SDjFJ0hhb1kuOSpJO3bJdclSSdHqW5ZKjkqTTt1yXHJUknaYmlxw1ICRpzG3eMD2SQDiRQ0ySpKJWH0E0NflEkk6GE+VGrMnJJ5I0KCfKNaDJySeSNCgnyjWgycknkjQoJ8o1oMnJJ5I0KCfKNaDJySeSNCgnyjWgycknkjSoJj+rIjNr30ldOp1OdrvdpsuQpIkSEfszs7PUdq0dYpIk9WdASJKKDAhJUpEBIUkqMiAkSUUGhCSpyICQJBUZEJKkIgNCklRkQEiSimoLiIi4NSKORsRDC9o+GBEzEfH16uvqBe/dHBGHIuJgRGysqy5J0mDqvFnfbcCfA7ef0P7RzPyzhQ0RcTFwHfDzwPnA30fEhZl5HElquaaWHK3tCCIz7wWeHHDzTcDnMvNHmfkd4BDwxrpqk6RJMb/k6MyxOZLnlxzdfWCm9n03cQ7ivRHxQDUEdVbVNg08tmCbw1WbJLVam5Yc/QTws8AbgMeBD1ftUdi2eB/yiNgWEd2I6M7OztZTpSSNidYsOZqZT2Tm8cz8MfApnh9GOgysXbDpBcCRRX7GzszsZGZnamqq3oIlqWGtWXI0Is5b8PJtwPwVTnuA6yLiZRHxamA98K+jrE2SxtGyXHI0Iu4ALgfWRMRh4APA5RHxBnrDR48AvwmQmQ9HxJ3AN4BngRu9gkmSXHL0lLnkqCSdPJcclSSdFgNCklRkQEiSigwISVKRASFJKjIgJElFBoQkqciAkCQVGRCSpCIDQpJUZEBIkooMCElSkQEhSSoyICRJRQaEJKnIgJAkFRkQkqQiA0KSVGRASJKKDAhJUpEBIUkqMiAkSUUGhCSpyICQJBUZEJKkokUDIiLuiYh1oytFkjRO+h1B3AZ8OSLeHxErR1SPJGlMnLHYG5l5Z0R8CfhDoBsRnwF+vOD9j4ygPklSQxYNiMozwP8ALwPOZEFASJKWt0UDIiKuBD4C7AEuycynR1aVJKlx/Y4g3g9cm5kPj6qYUdt9YIYdew9y5Ngc569exfaNF7F5w3TTZUnSCzT1WdXvHMQv1773Bu0+MMPNdz3I3DPHAZg5NsfNdz0IYEhIGhtNfla1dh7Ejr0Hn/sHnzf3zHF27D3YUEWS9GJNfla1NiCOHJs7qXZJakKTn1WtDYjzV686qXZJakKTn1WtDYjtGy9i1coVL2hbtXIF2zde1FBFkvRiTX5WLTUPYtmaP7njVUySxlmTn1WRmbXvpC6dTie73W7TZUjSRImI/ZnZWWq71g4xSZL6qy0gIuLWiDgaEQ8taDs7Ir4SEd+qHs+q2iMiPhYRhyLigYi4pK66JEmDqfMI4jbgyhPabgL2ZeZ6YF/1GuAqYH31tQ34RI11SZIGUFtAZOa9wJMnNG8CdlXPdwGbF7Tfnj33Aasj4ry6apMkLW3U5yBekZmPA1SP51bt08BjC7Y7XLW9SERsi4huRHRnZ2drLVaS2mxcTlJHoa14eVVm7szMTmZ2pqamai5Lktpr1AHxxPzQUfV4tGo/DKxdsN0FwJER1yZJWmDUAbEH2Fo93wrcvaD9+upqpkuBp+aHoiRJzahtJnVE3AFcDqyJiMPAB4BbgDsj4gbgUeDaavN7gKuBQ8DTwLvrqkuSNJjaAiIztyzy1hWFbRO4sa5aJEknb1xOUkuSxkxrb9YnSZNi7JYclSQ1zyVHJUlFLjkqSSpyyVFJUpFLjkqSilxyVJJU1OSSowaEJI25zRumRxIIJ3KISZJU1OojiKYmn0jSJGhtQDQ5+USSJkFrh5ianHwiSZOgtQHR5OQTSZoErQ2IJiefSNIkaG1ANDn5RJImQWtPUjc5+USSJkFrAwKam3wiSZOgtUNMkqT+Wn0E4UQ5SVpcawPCiXKS1F9rh5icKCdJ/bU2IJwoJ0n9tTYgnCgnSf21NiCcKCdJ/bX2JLUT5SSpv9YGBDhRTpL6ae0QkySpPwNCklRkQEiSigwISVKRASFJKjIgJElFBoQkqciAkCQVGRCSpCIDQpJUZEBIkooauRdTRDwC/BA4DjybmZ2IOBv4PLAOeAR4R2b+oM46XHJU0iRo6rOqySOIN2fmGzKzU72+CdiXmeuBfdXr2swvOTpzbI7k+SVHdx+YqXO3knRSmvysGqchpk3Arur5LmBznTtzyVFJk6DJz6qmAiKBL0fE/ojYVrW9IjMfB6gezy19Y0Rsi4huRHRnZ2dPuQCXHJU0CZr8rGoqIC7LzEuAq4AbI+JNg35jZu7MzE5mdqampk65AJcclTQJmvysaiQgMvNI9XgU+CLwRuCJiDgPoHo8WmcNLjkqaRI0+Vk18oCIiJ+KiDPnnwNvAR4C9gBbq822AnfXWcfmDdN86O2vZ3r1KgKYXr2KD7399V7FJGmsNPlZFZlZ+05esMOI19A7aoDeZbafzcw/iYhzgDuBVwKPAtdm5pP9flan08lut1trvZK03ETE/gVXkC5q5PMgMvPbwC8U2r8PXDHqeiRJZeN0maskaYwYEJKkIgNCklRkQEiSigwISVKRASFJKjIgJElFBoQkqciAkCQVGRCSpKJGlhwdFy45KkmLa21AzC/jN79S0/wyfoAhIUm0eIjJJUclqb/WBoRLjkpSf60NCJcclaT+WhsQLjkqSf219iT1/Ilor2KSpLLWBgT0QsJAkKSy1g4xSZL6a/URhBPlJGlxrQ0IJ8pJUn+tHWJyopwk9dfagHCinCT119qAcKKcJPXX2oBwopwk9dfak9ROlJOk/lobEOBEOUnqp7VDTJKk/gwISVJRq4eYnEktSYtrbUA4k1qS+mvtEJMzqSWpv9YGhDOpJam/1gaEM6klqb/WBoQzqSWpv9aepHYmtST119qAAGdSS1I/rR1ikiT1N3YBERFXRsTBiDgUETc1XY8ktdVYBURErAA+DlwFXAxsiYiLm61KktpprAICeCNwKDO/nZn/B3wO2NRwTZLUSuMWENPAYwteH67anhMR2yKiGxHd2dnZkRYnSW0yblcxRaEtX/AicyewEyAiZiPiu0PY7xrge0P4OZPC/i5vbepvm/oKw+vvqwbZaNwC4jCwdsHrC4Aji22cmVPD2GlEdDOzM4yfNQns7/LWpv62qa8w+v6O2xDTvwHrI+LVEfFS4DpgT8M1SVIrjdURRGY+GxHvBfYCK4BbM/PhhsuSpFYaq4AAyMx7gHtGvNudI95f0+zv8tam/raprzDi/kZmLr2VJKl1xu0chCRpTLQqIJa6jUdEvCwiPl+9f39ErBt9lcMzQH9/PyK+EREPRMS+iBjo0rdxNOgtWiLimojIiJjoK18G6W9EvKP6/T4cEZ8ddY3DNMDf8isj4msRcaD6e766iTqHISJujYijEfHQIu9HRHys+rd4ICIuqa2YzGzFF72T3v8FvAZ4KfAfwMUnbPNbwCer59cBn2+67pr7+2bgJ6vn75nU/g7S12q7M4F7gfuATtN11/y7XQ8cAM6qXp/bdN0193cn8J7q+cXAI03XfRr9fRNwCfDQIu9fDfwtvXljlwL311VLm44gBrmNxyZgV/X8r4ErIqI0eW8SLNnfzPxaZj5dvbyP3ryTSTToLVr+GPhT4H9HWVwNBunvbwAfz8wfAGTm0RHXOEyD9DeBn66e/wx95k+Nu8y8F3iyzyabgNuz5z5gdUScV0ctbQqIJW/jsXCbzHwWeAo4ZyTVDd8g/V3oBnr/VzKJBrlFywZgbWb+zSgLq8kgv9sLgQsj4p8j4r6IuHJk1Q3fIP39IPDOiDhM7yrI3x5NaY042f+2T9nYXeZaoyVv4zHgNpNi4L5ExDuBDvArtVZUn759jYiXAB8F3jWqgmo2yO/2DHrDTJfTOzL8x4h4XWYeq7m2OgzS3y3AbZn54Yj4ReAzVX9/XH95Izeyz6k2HUEMchuP57aJiDPoHar2O9QbZwPdtiQifhV4P/DWzPzRiGobtqX6eibwOuAfIuIReuO2eyb4RPWgf8t3Z+Yzmfkd4CC9wJhEg/T3BuBOgMz8F+An6N23aDk6qVsSnY42BcQgt/HYA2ytnl8DfDWrs0ITaMn+VsMuf0EvHCZ5jLpvXzPzqcxck5nrMnMdvfMtb83MbjPlnrZB/pZ307sIgYhYQ2/I6dsjrXJ4Bunvo8AVABHxWnoBsVxv97wHuL66mulS4KnMfLyOHbVmiCkXuY1HRPwR0M3MPcCn6R2aHqJ35HBdcxWfngH7uwN4OfBX1bn4RzPzrY0VfYoG7OuyMWB/9wJviYhvAMeB7Zn5/eaqPnUD9vd9wKci4vfoDbe8a1L/5y4i7qA3NLimOqfyAWAlQGZ+kt45lquBQ8DTwLtrq2VC/w0lSTVr0xCTJOkkGBCSpCIDQpJUZEBIkooMCElSkQEhDVlEvK26Y+zPNV2LdDoMCGn4tgD/xATPo5HAeRDSUEXEy+nd1uLNwJ7M9ChCE8sjCGm4NgN/l5n/CTxZ62IuUs0MCGm4ttBbr4DqcUuDtUinxSEmaUgi4hx6d9o8Su9+QCuqx1dN6n2B1G4eQUjDcw29lb5eVd05di3wHeCXGq5LOiUGhDQ8W4AvntD2BeDXGqhFOm0OMUmSijyCkCQVGRCSpCIDQpJUZEBIkooMCElSkQEhSSoyICRJRQaEJKno/wG/VMQKZu8K2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(df1['A'], df1['Y'])\n",
    "ax.set_xlabel('A')\n",
    "ax.set_ylabel('Y');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"8\" halign=\"left\">Y</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8.0</td>\n",
       "      <td>67.50</td>\n",
       "      <td>53.117121</td>\n",
       "      <td>10.0</td>\n",
       "      <td>27.5</td>\n",
       "      <td>60.0</td>\n",
       "      <td>87.5</td>\n",
       "      <td>170.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8.0</td>\n",
       "      <td>146.25</td>\n",
       "      <td>58.294205</td>\n",
       "      <td>50.0</td>\n",
       "      <td>105.0</td>\n",
       "      <td>160.0</td>\n",
       "      <td>185.0</td>\n",
       "      <td>220.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Y                                                     \n",
       "  count    mean        std   min    25%    50%    75%    max\n",
       "A                                                           \n",
       "0   8.0   67.50  53.117121  10.0   27.5   60.0   87.5  170.0\n",
       "1   8.0  146.25  58.294205  50.0  105.0  160.0  185.0  220.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df1.groupby('A').describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Now suppose treatment A is a polytomous variable that can take 4 possible values\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2 = pd.DataFrame({\n",
    "    'A': [\n",
    "        1, 1, 1, 1, 2, 2, 2, 2,\n",
    "        3, 3, 3, 3, 4, 4, 4, 4\n",
    "    ],\n",
    "    'Y': [\n",
    "        110, 80, 50, 40, 170, 30, 70, 50,\n",
    "        110, 50, 180, 130, 200, 150, 220, 210\n",
    "    ]\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFoBJREFUeJzt3X+QXWd93/H3p0Kl28JUgBdiyxYCxqj8CJXNjscdF+KEtDIOxYZCarcFQ+gItzANU0YFQSZQ/glTFUhpZmBEYDApcexiIVxiIhxISpipSVfIsUyNwCYGtNJYwkKY1juuJX/7xz1rX6/PalfW3j337n2/Znb23uece8/38eOdj86P55xUFZIkzfc3ui5AkjScDAhJUisDQpLUyoCQJLUyICRJrQwISVIrA0KS1MqAkCS1MiAkSa2e0nUBZ+Kss86qjRs3dl2GJI2UvXv3/qSqJhdbb2ABkeQ84HPALwCPADur6j8n2QH8E+D/AfcAb62q40k2AncBB5qvuK2qrj3VNjZu3Mj09PSAeiBJq1OSHy5lvUEeYjoBvLuqXgRcDLwjyYuBW4GXVtXLgO8B2/s+c09VbW5+ThkOkqTBGlhAVNXhqvp28/rn9PYO1lfVV6vqRLPabcC5g6pBkvTkrchJ6ubw0QXAt+Yt+g3gK33vn5dkX5L/keQVK1GbJKndwE9SJ3kacBPwrqp6oK/9/fQOQ32+aToMbKiq+5O8HNid5CX9n2k+txXYCrBhw4ZBly9JY2ugexBJ1tILh89X1a6+9muA1wD/opoHUlTVQ1V1f/N6L70T2C+c/51VtbOqpqpqanJy0ZPwkqQnaZBXMQX4NHBXVX20r/0y4D3AL1XVg33tk8CxqjqZ5PnA+cAPBlWfJI2a3ftm2LHnAIeOz3LOugm2bdnElResH9j2BnmI6RLgTcD+JLc3be8DPg48Fbi1lyGPXs76SuBDSU4AJ4Frq+rYAOuTpJGxe98M23ftZ/bhkwDMHJ9l+679AAMLiYEFRFV9E0jLolsWWP8meoejJEnz7Nhz4NFwmDP78El27DkwsIDwVhuSNAIOHZ89rfblYEBI0gg4Z93EabUvBwNCkkbAti2bmFi75nFtE2vXsG3LpoFtc6Rv1idJ42LuPMNquYpJkrSMrrxg/UADYT4PMUmSWhkQkqRWBoQkqZXnICRpRKymW21IkpZJF7fa8BCTJI2AU91qY1AMCEkaAd5qQ5LUylttSJJaeasNSVIrb7UhSVqQt9qQJA2FgQVEkvOS/FmSu5J8J8lvNu3PTHJrku83v5/RtCfJx5PcneSOJBcOqjZJ0uIGuQdxAnh3Vb0IuBh4R5IXA+8FvlZV5wNfa94DvBo4v/nZCnxigLVJ0sjZvW+GSz78dZ733j/mkg9/nd37Zga6vYEFRFUdrqpvN69/DtwFrAeuAK5rVrsOuLJ5fQXwueq5DViX5OxB1SdJo2RuJvXM8VmKx2ZSDzIkVuQcRJKNwAXAt4DnVNVh6IUI8OxmtfXAj/s+drBpk6SxtypnUid5GnAT8K6qeuBUq7a0Vcv3bU0ynWT66NGjy1WmJA21VTeTOslaeuHw+ara1TTfN3foqPl9pGk/CJzX9/FzgUPzv7OqdlbVVFVNTU5ODq54SRoiq2omdZIAnwbuqqqP9i26GbimeX0N8KW+9jc3VzNdDPxs7lCUJI271TaT+hLgTcD+JLc3be8DPgzcmORtwI+ANzbLbgEuB+4GHgTeOsDaJGmkdDGTOlVPOMw/Mqampmp6errrMiRppCTZW1VTi63nTGpJUisDQpLUyoCQJLUyICRJrQwISVIrA0KS1MqAkCS1MiAkSa0MCElSKwNCktTKgJAktTIgJEmtDAhJUqtB3u5b0gjbvW9mRW8treFjQEh6gt37Zti+a/+jz0CeOT7L9l37AQyJMeIhJklPsGPPgUfDYc7swyfZsedARxWpCwaEpCc4dHz2tNq1Og3ymdSfSXIkyZ19bTckub35uXfuUaRJNiaZ7Vv2yUHVJWlx56ybOK12rU6D3IP4LHBZf0NV/bOq2lxVm4GbgF19i++ZW1ZV1w6wLkmL2LZlExNr1zyubWLtGrZt2dRRRerCwE5SV9U3kmxsW5YkwK8DvzKo7Ut68uZORHsV03jr6iqmVwD3VdX3+9qel2Qf8ADwW1X1F92UJgl6IWEgjLeuAuJq4Pq+94eBDVV1f5KXA7uTvKSqHpj/wSRbga0AGzZsWJFiJWkcrfhVTEmeArweuGGuraoeqqr7m9d7gXuAF7Z9vqp2VtVUVU1NTk6uRMmSNJa62IP4VeC7VXVwriHJJHCsqk4meT5wPvCDDmpTR5y1Kw2fQV7mej3wP4FNSQ4meVuz6Coef3gJ4JXAHUn+CvgCcG1VHRtUbRouc7N2Z47PUjw2a3f3vpmuS5PG2iCvYrp6gfa3tLTdRO+yV42hU83adS9C6o4zqdU5Z+1Kw8mAUOectSsNJwNCnXPWrjScvN23OuesXWk4GRAaCs7alYaPh5gkSa0MCElSKwNCktTKgJAktTIgJEmtDAhJUisDQpLUyoCQJLUyICRJrQwISVIrA0KS1MqAkCS1GuQjRz+T5EiSO/vaPphkJsntzc/lfcu2J7k7yYEkWwZVlyRpaQZ5N9fPAr8HfG5e+8eq6j/1NyR5Mb1nVb8EOAf40yQvrKqTSJKA3vPbV/K2+APbg6iqbwDHlrj6FcAfVdVDVfXXwN3ARYOqTZJGze59M2zftZ+Z47MUMHN8lu279rN738zAttnFOYh3JrmjOQT1jKZtPfDjvnUONm2SJHoP1Jp9+PEHVWYfPsmOPQcGts2VDohPAC8ANgOHgY807WlZt9q+IMnWJNNJpo8ePTqYKiVpyBw6Pnta7cthRQOiqu6rqpNV9QjwKR47jHQQOK9v1XOBQwt8x86qmqqqqcnJycEWLElD4px1E6fVvhxWNCCSnN339nXA3BVONwNXJXlqkucB5wN/uZK1SdIw27ZlExNr1zyubWLtGrZt2TSwbQ7sKqYk1wOXAmclOQh8ALg0yWZ6h4/uBd4OUFXfSXIj8L+BE8A7vIJJkh4zd7XSSl7FlKrWQ/0jYWpqqqanp7suQ5JGSpK9VTW12HrOpJYktTIgJEmtDAhJUisDQpLUyoCQJLUyICRJrQwISVIrA0KS1MqAkCS1MiAkSa0MCElSKwNCktTKgJAktRrY7b4ljbbd+2ZW9NbSGj4GhKQn2L1vhu279j/6DOSZ47Ns37UfwJAYIx5ikvQEO/YceDQc5sw+fJIdew50VJG6YEBIeoJDx2dPq12r08ACIslnkhxJcmdf244k301yR5IvJlnXtG9MMpvk9ubnk4OqS9Lizlk3cVrtWp0WDIgktyTZeAbf/VngsnlttwIvraqXAd8Dtvctu6eqNjc/157BdiWdoW1bNjGxds3j2ibWrmHblk0dVaQunGoP4rPAV5O8P8na0/3iqvoGcGxe21er6kTz9jbg3NP9XkmDd+UF6/md1/8i69dNEGD9ugl+5/W/6AnqMbPgVUxVdWOSPwZ+G5hO8gfAI33LP3qG2/4N4Ia+989Lsg94APitqvqLM/x+SWfgygvWGwhjbrHLXB8G/i/wVODp9AXEmUjyfuAE8Pmm6TCwoaruT/JyYHeSl1TVAy2f3QpsBdiwYcNylCNJarFgQCS5DPgocDNwYVU9uBwbTHIN8BrgVVVVAFX1EPBQ83pvknuAFwLT8z9fVTuBnQBTU1O1HDVJkp7oVHsQ7wfeWFXfWa6NNaHzHuCX+gMnySRwrKpOJnk+cD7wg+XariTp9J3qHMQrzuSLk1wPXAqcleQg8AF6Vy09Fbg1CcBtzRVLrwQ+lOQEcBK4tqqOtX6xJGlFDOxWG1V1dUvzpxdY9ybgpkHVIkk6fc6kliS1MiAkSa0MCElSKwNCktTKgJAktTIgJEmtDAhJUisDQpLUaiyfSe3D2KXF+XeisQsIH8YuLc6/E8EYHmLyYezS4vw7EYxhQPgwdmlx/p0IxjAgfBi7tDj/TgRjGBA+jF1anH8ngjE8ST13gs2rM6SF+XcigDRP/RxJU1NTNT39hKeSSpJOIcneqppabL2xO8QkSVqagQZEks8kOZLkzr62Zya5Ncn3m9/PaNqT5ONJ7k5yR5ILB1mbJOnUBr0H8Vngsnlt7wW+VlXnA19r3gO8Gji/+dkKfGLAtUmSTmGgAVFV3wCOzWu+AriueX0dcGVf++eq5zZgXZKzB1mfJGlhXZyDeE5VHQZofj+7aV8P/LhvvYNNmySpA8N0kjotbU+4xCrJ1iTTSaaPHj26AmVJ0njqIiDumzt01Pw+0rQfBM7rW+9c4ND8D1fVzqqaqqqpycnJgRcrSeOqi4C4GbimeX0N8KW+9jc3VzNdDPxs7lCUJGnlDXQmdZLrgUuBs5IcBD4AfBi4McnbgB8Bb2xWvwW4HLgbeBB46yBrkySd2kADoqquXmDRq1rWLeAdg6xHkrR0w3SSWpI0RAwISVIrA0KS1MqAkCS1MiAkSa0MCElSKwNCktTKgJAktTIgJEmtDAhJUisDQpLUyoCQJLUyICRJrQZ6N9dhtXvfDDv2HODQ8VnOWTfBti2buPICn24qSf3GLiB275th+679zD58EoCZ47Ns37UfwJCQpD5jd4hpx54Dj4bDnNmHT7Jjz4GOKpKk4TR2AXHo+OxptUvSuFrxgEiyKcntfT8PJHlXkg8mmelrv3wQ2z9n3cRptUvSuFrxgKiqA1W1uao2Ay+n9/zpLzaLPza3rKpuGcT2t23ZxMTaNY9rm1i7hm1bNg1ic5I0sro+Sf0q4J6q+mGSFdng3Ilor2KSpFPrOiCuAq7ve//OJG8GpoF3V9VPB7HRKy9YbyBI0iI6O0md5G8CrwX+W9P0CeAFwGbgMPCRBT63Ncl0kumjR4+uSK2SNI66vIrp1cC3q+o+gKq6r6pOVtUjwKeAi9o+VFU7q2qqqqYmJydXsFxJGi9dHmK6mr7DS0nOrqrDzdvXAXd2UpU64ex2afh0EhBJ/jbwj4C39zX/xySbgQLunbdMq5iz26Xh1ElAVNWDwLPmtb2pi1rUvVPNbjcgpO6M3UxqDR9nt0vDyYBQ55zdLg0nA0Kdc3a7NJy6nignObtdGlIGhIaCs9ul4eMhJklSKwNCktTKgJAktTIgJEmtDAhJUisDQpLUyoCQJLUyICRJrQwISVIrA0KS1MqAkCS1MiAkSa06u1lfknuBnwMngRNVNZXkmcANwEZ6jx399ar6aVc1StI46/purr9cVT/pe/9e4GtV9eEk723ev2e5N7p734y3lh4yjsnwcUzUdUDMdwVwafP6OuDPWeaA2L1vhu279j/6DOSZ47Ns37UfwP/5O+KYDB/HRNDtOYgCvppkb5KtTdtzquowQPP72cu90R17Djz6P/2c2YdPsmPPgeXelJbIMRk+jomg2z2IS6rqUJJnA7cm+e5SPtSEyVaADRs2nPZGDx2fPa12DZ5jMnwcE0GHexBVdaj5fQT4InARcF+SswGa30daPrezqqaqampycvK0t3vOuonTatfgOSbDxzERdBQQSf5OkqfPvQb+MXAncDNwTbPaNcCXlnvb27ZsYmLtmse1Taxdw7Ytm5Z7U1oix2T4OCaC7g4xPQf4YpK5Gv6wqv4kyf8CbkzyNuBHwBuXe8NzJ9i8OmN4OCbDxzERQKqq6xqetKmpqZqenu66DEkaKUn2VtXUYus5k1qS1MqAkCS1GraJcivCGaKStLixCwhniErS0ozdISZniErS0oxdQDhDVJKWZuwCwhmikrQ0YxcQzhCVpKUZu5PUzhCVpKUZu4CAXkgYCJJ0amN3iEmStDQGhCSplQEhSWo1lucgNHy8/Yk0fAwIdc7bn0jDyUNM6py3P5GGkwGhznn7E2k4rXhAJDkvyZ8luSvJd5L8ZtP+wSQzSW5vfi5f6drUDW9/Ig2nLvYgTgDvrqoXARcD70jy4mbZx6pqc/NzSwe1qQPe/kQaTit+krqqDgOHm9c/T3IX4JnIMebtT6Th1OlVTEk2AhcA3wIuAd6Z5M3ANL29jJ92V51Wkrc/kYZPZyepkzwNuAl4V1U9AHwCeAGwmd4exkcW+NzWJNNJpo8ePbpi9UrSuOkkIJKspRcOn6+qXQBVdV9VnayqR4BPARe1fbaqdlbVVFVNTU5OrlzRkjRmuriKKcCngbuq6qN97Wf3rfY64M6Vrk2S9JguzkFcArwJ2J/k9qbtfcDVSTYDBdwLvL2D2iRJjS6uYvomkJZFXtYqSUMkVdV1DU9akqPAD8/gK84CfrJM5XRttfRltfQDVk9fVks/YPX05Uz78dyqWvQk7kgHxJlKMl1VU13XsRxWS19WSz9g9fRltfQDVk9fVqof3otJktTKgJAktRr3gNjZdQHLaLX0ZbX0A1ZPX1ZLP2D19GVF+jHW5yAkSQsb9z0ISdICxiIgknwmyZEkrbOz0/PxJHcnuSPJhStd41IsoR+XJvlZ3zM1fnula1yKhZ4JMm+doR+TJfZjVMbkbyX5yyR/1fTlP7Ss89QkNzRj8q3mZptDZYn9eEuSo31j8q+6qHWpkqxJsi/Jl1uWDXZMqmrV/wCvBC4E7lxg+eXAV+hN4LsY+FbXNT/JflwKfLnrOpfQj7OBC5vXTwe+B7x41MZkif0YlTEJ8LTm9Vp6d1i+eN46/wb4ZPP6KuCGrut+kv14C/B7Xdd6Gn36d8Aftv1/NOgxGYs9iKr6BnDsFKtcAXyuem4D1s27N9RQWEI/RkJVHa6qbzevfw60PRNk6Mdkif0YCc1/5//TvF3b/Mw/QXkFcF3z+gvAq5p7qw2NJfZjZCQ5F/g14PcXWGWgYzIWAbEE64Ef970/yIj+oQP/oNm9/kqSl3RdzGLmPROk30iNySn6ASMyJs2hjNuBI8CtVbXgmFTVCeBnwLNWtsrFLaEfAP+0OXT5hSTnrXCJp+N3gX8PPLLA8oGOiQHR05a4o/ivjm/Tm0L/94H/AuzuuJ5TankmyOMWt3xkKMdkkX6MzJhU73b7m4FzgYuSvHTeKiMxJkvox38HNlbVy4A/5bF/gQ+VJK8BjlTV3lOt1tK2bGNiQPQcBPr/FXEucKijWp60qnpgbve6es/0XpvkrI7LatX2TJB5RmJMFuvHKI3JnKo6Dvw5cNm8RY+OSZKnAH+XIT7kuVA/qur+qnqoefsp4OUrXNpSXQK8Nsm9wB8Bv5Lkv85bZ6BjYkD03Ay8ubly5mLgZ9V7dvZISfILc8cfk1xEb3zv77aqJ1romSDzDP2YLKUfIzQmk0nWNa8ngF8FvjtvtZuBa5rXbwC+Xs3Z0WGxlH7MO5f1WnrnjoZOVW2vqnOraiO9E9Bfr6p/OW+1gY5Jp8+kXilJrqd3NclZSQ4CH6B38oqq+iS9W41fDtwNPAi8tZtKT20J/XgD8K+TnABmgauG7Q+4sdAzQTbASI3JUvoxKmNyNnBdkjX0QuzGqvpykg8B01V1M70w/IMkd9P7V+pV3ZW7oKX0498meS1wgl4/3tJZtU/CSo6JM6klSa08xCRJamVASJJaGRCSpFYGhCSplQEhSWplQEjLLMnrklSSv9d1LdKZMCCk5Xc18E2Gc56AtGTOg5CWUXNfpgPALwM3V5V7ERpZ7kFIy+tK4E+q6nvAsWF80JG0VAaEtLyupndjNZrfV3dYi3RGPMQkLZMkz6J3d80j9G65vKb5/dwhvf+SdEruQUjL5w30noL33KraWFXnAX8N/MOO65KeFANCWj5XA1+c13YT8M87qEU6Yx5ikiS1cg9CktTKgJAktTIgJEmtDAhJUisDQpLUyoCQJLUyICRJrQwISVKr/w8TcQgpoTjBBwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(df2['A'], df2['Y'])\n",
    "ax.set_xlabel('A')\n",
    "ax.set_ylabel('Y');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"8\" halign=\"left\">Y</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>31.622777</td>\n",
       "      <td>40.0</td>\n",
       "      <td>47.5</td>\n",
       "      <td>65.0</td>\n",
       "      <td>87.5</td>\n",
       "      <td>110.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.0</td>\n",
       "      <td>80.0</td>\n",
       "      <td>62.182527</td>\n",
       "      <td>30.0</td>\n",
       "      <td>45.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>95.0</td>\n",
       "      <td>170.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.0</td>\n",
       "      <td>117.5</td>\n",
       "      <td>53.774219</td>\n",
       "      <td>50.0</td>\n",
       "      <td>95.0</td>\n",
       "      <td>120.0</td>\n",
       "      <td>142.5</td>\n",
       "      <td>180.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.0</td>\n",
       "      <td>195.0</td>\n",
       "      <td>31.091264</td>\n",
       "      <td>150.0</td>\n",
       "      <td>187.5</td>\n",
       "      <td>205.0</td>\n",
       "      <td>212.5</td>\n",
       "      <td>220.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Y                                                     \n",
       "  count   mean        std    min    25%    50%    75%    max\n",
       "A                                                           \n",
       "1   4.0   70.0  31.622777   40.0   47.5   65.0   87.5  110.0\n",
       "2   4.0   80.0  62.182527   30.0   45.0   60.0   95.0  170.0\n",
       "3   4.0  117.5  53.774219   50.0   95.0  120.0  142.5  180.0\n",
       "4   4.0  195.0  31.091264  150.0  187.5  205.0  212.5  220.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df2.groupby('A').describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 11.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 11.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "A, Y = zip(*(\n",
    "    (3, 21),\n",
    "    (11, 54),\n",
    "    (17, 33),\n",
    "    (23, 101),\n",
    "    (29, 85),\n",
    "    (37, 65),\n",
    "    (41, 157),\n",
    "    (53, 120),\n",
    "    (67, 111),\n",
    "    (79, 200),\n",
    "    (83, 140),\n",
    "    (97, 220),\n",
    "    (60, 230),\n",
    "    (71, 217),\n",
    "    (15, 11),\n",
    "    (45, 190),\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "df3 = pd.DataFrame({'A': A, 'Y': Y, 'constant': np.ones(16)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEbpJREFUeJzt3X+MHGd9x/H3t06gx4/qEnKJ4ouDQ+Ua0lrE0QmlDa0CqeqEovqwSJu0FVZE5f4RVKhoKodKpf2jctrQ0CKqSGlJ41Q0EIFxrBJhUgcJ+IOUM67yg+BiQZr47MamwYESCxzn2z92lmzsG3vPt7OzO/N+SafdfW68993duft4npnneSIzkSRpIT9TdwGSpNFlSEiSShkSkqRShoQkqZQhIUkqZUhIkkoZEpKkUoaEJKmUISFJKnVW3QUsxXnnnZcrV66suwxJGiu7d+/+XmZO9bPtWIfEypUrmZubq7sMSRorEfHf/W5rd5MkqZQhIUkqZUhIkkoZEpKkUoaEJKnUWF/dJOlk2/fMc9vOvRw4cpTlkxPcvG41s2un6y5LY8qQkBpk+555btn2KEePHQdg/shRbtn2KIBBoTNid5PUILft3PvTgOg6euw4t+3cW1NFGneGhNQgB44cXVS7dDqGhNQgyycnFtUunY4hITXIzetWM3H2spe1TZy9jJvXra6pIo07T1xLDdI9Oe3VTRoUQ0JqmNm104aCBsbuJklSKY8kJGmMDHuwpCEhSWOijsGSdjdJ0pioY7CkRxKSKuVcUoNTx2BJjyQkVabbPTJ/5CjJS90j2/fM113aWKpjsKQhIakyziU1WHUMlrS7SVJlnEtqsOoYLGlISKrM8skJ5hcIBOeSOnPDHixpd5OkyjiX1PjzSEJSZZxLavwZEpIq5VxS483uJklSKUNCklTKkJAklTIkJEmlDAlJUilDQpJUypCQJJWqLCQiYkVEfCkinoiIxyPi/UX7uRHxYER8u7g9p2iPiPhYROyLiEci4vKqapPULNv3zHPlrQ9xyebPc+WtDznL7ABVeSTxAvDBzHwTcAVwU0RcCmwGdmXmKmBX8RjgWmBV8bUJuKPC2iQ1hNORV6uykMjMg5n5jeL+D4EngGlgPbC12GwrMFvcXw/ckx1fAyYj4sKq6pPUDE5HXq2hnJOIiJXAWuBh4ILMPAidIAHOLzabBp7u+Wf7i7YTn2tTRMxFxNzhw4erLFvSGHA68mpVHhIR8Rrgs8AHMvMHp9p0gbY8qSHzzsycycyZqampQZUpaUzVsVpbm1QaEhFxNp2A+GRmbiuan+l2IxW3h4r2/cCKnn9+EXCgyvokjT+nI69WlVc3BfAJ4InMvL3nWzuAjcX9jcD9Pe3vKa5yugJ4rtstJUllZtdOs2XDGqYnJwhgenKCLRvWOPPsgETmST06g3niiLcCXwEeBV4smj9E57zEfcDFwFPAdZn5bBEqHweuAZ4HbszMuVP9jJmZmZybO+UmkqQTRMTuzJzpZ9vK1pPIzK+y8HkGgKsX2D6Bm6qqR5K0eC46JPXYvmfeVdSkHoaEVOgOyupec98dlAUYFGot526SCg7Kkk5mSEgFB2VJJzMkpIKDsqSTGRJSwUFZ0sk8cS0VuienvbpJeokhIfWYXTttKEg97G6SJJUyJCRJpQwJSVIpQ0KSVMqQkCSVMiQkSaUMCUlSKUNCklTKkJAklTIkJEmlDAlJUilDQpJUypCQJJUyJCRJpQwJSVIpQ0KSVMqQkCSVMiQkSaUMCUlSKUNCklTKkJAklTIkJEmlzqq7AKkq2/fMc9vOvRw4cpTlkxPcvG41s2un6y5LGiuGhBpp+555btn2KEePHQdg/shRbtn2KIBBIS2C3U1qpNt27v1pQHQdPXac23burakiaTwZEmqkA0eOLqpd0sIMCTXS8smJRbVLWpghoUa6ed1qJs5e9rK2ibOXcfO61TVVJI2nykIiIu6KiEMR8VhP219ExHxE/Gfx9Y6e790SEfsiYm9ErKuqLrXD7NpptmxYw/TkBAFMT06wZcMaT1pLi1Tl1U13Ax8H7jmh/aOZ+ZHehoi4FLge+EVgOfDvEfELmXkc6QzNrp02FKQlquxIIjO/DDzb5+brgU9l5o8z87vAPuAtVdUmSepPHeck3hcRjxTdUecUbdPA0z3b7C/aThIRmyJiLiLmDh8+XHWtktRqww6JO4CfBy4DDgJ/W7THAtvmQk+QmXdm5kxmzkxNTVVTpSQJGHJIZOYzmXk8M18E/pGXupT2Ayt6Nr0IODDM2iRJJxtqSETEhT0P3wV0r3zaAVwfEa+MiEuAVcB/DLM2SdLJKru6KSLuBa4CzouI/cCHgasi4jI6XUlPAn8IkJmPR8R9wDeBF4CbvLJJkuoXmQt2/Y+FmZmZnJubq7sMSRorEbE7M2f62dZZYCW1klPJ98eQkNQ6TiXfP+duktQ6TiXfP0NCUus4lXz/DAlJreNU8v0zJCS1jlPJ988T15Jap3ty2qubTs+QkNRKTiXfH7ubJEmlDAlJUilDQpJUypCQJJUyJCRJpQwJSVIpQ0KSVMqQkCSVMiQkSaUMCUlSqdKQiIgHImLl8EqRJI2aU83ddDfwxYjYCvxNZh4bTklS87l0psZFaUhk5n0R8Xngz4G5iPgX4MWe798+hPqkxnHpTI2T052TOAb8CHgl8NoTviSdAZfO1DgpPZKIiGuA24EdwOWZ+fzQqpIazKUzNU5OdU7iz4DrMvPxYRUjtcHyyQnmFwgEl87UKCrtbsrMXzUgpMFz6UyNE1emk4bMpTM1TgwJqQYunalxYUhI6pvjO9rHkJDUF8d3tJNzN0nqi+M72smQkNQXx3e0kyEhqS9l4zgc39FshoSkvji+o508cS2pL47vaCdDQlLfHN/RPnY3SZJKVRYSEXFXRByKiMd62s6NiAcj4tvF7TlFe0TExyJiX0Q8EhGXV1WXOte7X3nrQ1yy+fNceetDbN8zX3dJkkZUlUcSdwPXnNC2GdiVmauAXcVjgGuBVcXXJuCOCutqte6AqPkjR0leGhBlUEhaSGUhkZlfBp49oXk9sLW4vxWY7Wm/Jzu+BkxGxIVV1dZmDoiStBjDPidxQWYeBChuzy/ap4Gne7bbX7RpwBwQJWkxRuXEdSzQlgtuGLEpIuYiYu7w4cMVl9U8DoiStBjDDolnut1Ixe2hon0/sKJnu4uAAws9QWbemZkzmTkzNTVVabFN5IAoSYsx7JDYAWws7m8E7u9pf09xldMVwHPdbikN1uzaabZsWMP05AQBTE9OsGXDGq99l7SgygbTRcS9wFXAeRGxH/gwcCtwX0S8F3gKuK7Y/AHgHcA+4HngxqrqkgOiJPWvspDIzBtKvnX1AtsmcFNVtUiSzsyonLiWJI0gQ0KSVMqQkCSVMiQkSaUMCUlSKdeT0JJt3zPvQjRSQxkSWpLurLLdSQO7s8oCBoXUAHY3aUmcVVZqNkNCS+KsslKzGRJaEmeVlZrNkNCSOKus1GyeuNaSdE9Oe3WT1EyGhJbMWWWl5rK7SZJUypCQJJUyJCRJpQwJSVIpQ0KSVMqQkCSVMiQkSaUMCUlSKUNCklTKkJAklXJaDo0NV8CThs+Q0FhwBTypHnY3aSy4Ap5UD0NCY8EV8KR6GBIaC66AJ9XDkNBYcAU8qR6euNZYcAU8qR6GhMaGK+BJw2d3kySplCEhSSpld1NFHB0sqQkMiQo4OlhSU9jdVAFHB0tqilqOJCLiSeCHwHHghcyciYhzgU8DK4Engd/OzO/XUd9SOTpYUlPUeSTxtsy8LDNnisebgV2ZuQrYVTweS44OltQUo9TdtB7YWtzfCszWWMuSODpYUlPUFRIJfDEidkfEpqLtgsw8CFDcnl9TbUs2u3aaLRvWMD05QQDTkxNs2bDGk9aSxk5dVzddmZkHIuJ84MGI+Fa//7AIlU0AF198cVX1LZmjgyU1QS1HEpl5oLg9BHwOeAvwTERcCFDcHir5t3dm5kxmzkxNTQ2rZElqpaGHRES8OiJe270P/AbwGLAD2FhsthG4f9i1SZJero7upguAz0VE9+f/a2Z+ISK+DtwXEe8FngKuq6E2SVKPoYdEZn4HePMC7f8LXD3seiRJ5UbpElhJ0ogxJCRJpQwJSVIpQ0KSVMqQkCSVMiQkSaUMCUlSKVemGxMuhyqpDobEGHA5VEl1sbtpDLgcqqS6GBJjwOVQJdXFkBgDLocqqS6GxBhwOVRJdfHE9Rjonpz26iZJw2ZIjAmXQ5VUB7ubJEmlDAlJUilDQpJUqnXnJJzeQpL616qQcHoLSVqcVnU3Ob2FJC1Oq0LC6S0kaXFaFRJObyFJi9OqkHB6C0lanFaduHZ6C0lanFaFBDi9hSQtRqu6myRJi9O6I4lR5kA/SaPGkBgRDvSTNIrsbhoRDvSTNIoMiRHhQD9Jo8iQGBEO9JM0igyJEeFAP0mjyBPXI8KBfpJGkSExQhzoJ2nU2N0kSSo1ciEREddExN6I2BcRm+uuR5LabKRCIiKWAf8AXAtcCtwQEZfWW5UktddIhQTwFmBfZn4nM38CfApYX3NNktRaoxYS08DTPY/3F22SpBqM2tVNsUBbvmyDiE3ApuLh/0XEQvNWnAd8b8C1jZO2v37wPQDfA/A9KHv9r+/3CUYtJPYDK3oeXwQc6N0gM+8E7jzVk0TEXGbODL688dD21w++B+B7AL4Hg3j9o9bd9HVgVURcEhGvAK4HdtRckyS11kgdSWTmCxHxPmAnsAy4KzMfr7ksSWqtkQoJgMx8AHhgiU9zyu6oFmj76wffA/A9AN+DJb/+yMzTbyVJaqVROychSRohjQqJNk7pERErIuJLEfFERDweEe8v2s+NiAcj4tvF7Tl111qliFgWEXsi4t+Kx5dExMPF6/90cSFEY0XEZER8JiK+VewLv9zCfeCPi9+BxyLi3oj42abvBxFxV0QciojHetoW/Nyj42PF38dHIuLyfn5GY0KixVN6vAB8MDPfBFwB3FS87s3ArsxcBewqHjfZ+4Eneh7/NfDR4vV/H3hvLVUNz98DX8jMNwJvpvNetGYfiIhp4I+Amcz8JToXvlxP8/eDu4FrTmgr+9yvBVYVX5uAO/r5AY0JCVo6pUdmHszMbxT3f0jnj8M0nde+tdhsKzBbT4XVi4iLgN8E/ql4HMDbgc8UmzT99f8c8GvAJwAy8yeZeYQW7QOFs4CJiDgLeBVwkIbvB5n5ZeDZE5rLPvf1wD3Z8TVgMiIuPN3PaFJItH5Kj4hYCawFHgYuyMyD0AkS4Pz6Kqvc3wF/CrxYPH4dcCQzXygeN31feANwGPjnosvtnyLi1bRoH8jMeeAjwFN0wuE5YDft2g+6yj73M/ob2aSQOO2UHk0WEa8BPgt8IDN/UHc9wxIR7wQOZebu3uYFNm3yvnAWcDlwR2auBX5Eg7uWFlL0u68HLgGWA6+m071yoibvB6dzRr8XTQqJ007p0VQRcTadgPhkZm4rmp/pHkoWt4fqqq9iVwK/FRFP0ulifDudI4vJotsBmr8v7Af2Z+bDxePP0AmNtuwDAL8OfDczD2fmMWAb8Cu0az/oKvvcz+hvZJNCopVTehT9758AnsjM23u+tQPYWNzfCNw/7NqGITNvycyLMnMlnc/8ocz8PeBLwLuLzRr7+gEy83+ApyNiddF0NfBNWrIPFJ4CroiIVxW/E933oDX7QY+yz30H8J7iKqcrgOe63VKn0qjBdBHxDjr/i+xO6fFXNZdUuYh4K/AV4FFe6pP/EJ3zEvcBF9P5BbouM088wdUoEXEV8CeZ+c6IeAOdI4tzgT3A72fmj+usr0oRcRmdE/evAL4D3EjnP4Gt2Qci4i+B36Fzxd8e4A/o9Lk3dj+IiHuBq+jM9voM8GFgOwt87kV4fpzO1VDPAzdm5txpf0aTQkKSNFhN6m6SJA2YISFJKmVISJJKGRKSpFKGhCSplCEhDVhEvCsiMiLeWHct0lIZEtLg3QB8lc7gPmmsOU5CGqBiDq29wNuAHcXU3dLY8khCGqxZOus6/BfwbL8Lu0ijypCQBusGOtNAUNzeUGMt0pLZ3SQNSES8js5Mm4foTMG8rLh9ffqLpjHlkYQ0OO+ms/LX6zNzZWauAL4LvLXmuqQzZkhIg3MD8LkT2j4L/G4NtUgDYXeTJKmURxKSpFKGhCSplCEhSSplSEiSShkSkqRShoQkqZQhIUkqZUhIkkr9PwCMrhIE0jZ6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(df3.A, df3.Y)\n",
    "ax.set_xlabel('A')\n",
    "ax.set_ylabel('Y');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols = sm.OLS(Y, df3[['constant', 'A']])\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/pettychen/anaconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>   24.5464</td> <td>   21.330</td> <td>    1.151</td> <td> 0.269</td> <td>  -21.202</td> <td>   70.295</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>        <td>    2.1372</td> <td>    0.400</td> <td>    5.347</td> <td> 0.000</td> <td>    1.280</td> <td>    2.994</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = res.summary()\n",
    "summary.tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I couldn't figure out how to get confidence intervals for expected values from Statsmodels. I'll calculate them from scratch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = df3.shape[0]\n",
    "yvar = (res.resid * res.resid).sum() / (n - 2)  # = res.mse_resid\n",
    "xval = np.array([[1, 90]])\n",
    "X = df3[['constant', 'A']]\n",
    "XpXinv = np.linalg.inv(np.dot(X.T, X))\n",
    "se_mean = np.sqrt(yvar * np.dot(xval, np.dot(XpXinv, xval.T)))[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate      95% C.I.\n",
      "E[Y|A=90]   216.89   (172.15, 261.63)\n"
     ]
    }
   ],
   "source": [
    "t = scipy.stats.t.ppf(0.975, n - 2)\n",
    "ypred = res.predict([[1, 90]])[0]\n",
    "print('           estimate      95% C.I.')\n",
    "print(\n",
    "    'E[Y|A=90]   {:>6.2f}   ({:>6.2f}, {:>6.2f})'.format(\n",
    "        ypred, ypred - t * se_mean, ypred + t * se_mean\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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>mean</th>\n",
       "      <th>mean_se</th>\n",
       "      <th>mean_ci_lower</th>\n",
       "      <th>mean_ci_upper</th>\n",
       "      <th>obs_ci_lower</th>\n",
       "      <th>obs_ci_upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>216.890047</td>\n",
       "      <td>20.861403</td>\n",
       "      <td>172.146787</td>\n",
       "      <td>261.633308</td>\n",
       "      <td>112.271343</td>\n",
       "      <td>321.508751</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mean    mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  \\\n",
       "0  216.890047  20.861403     172.146787     261.633308    112.271343   \n",
       "\n",
       "   obs_ci_upper  \n",
       "0    321.508751  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Alternative method of calculating the expected value of Y when A =90\n",
    "# .get_prediction is an easier way of calculating the expected value and 95% CI\n",
    "pred = res.get_prediction(exog=[1,90])   \n",
    "pred.summary_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHZxJREFUeJzt3Xt4VOW59/HvzUGJVkmRF4QARqkbRCuiURF1C0iLIlas575V8UTd2m1rbRQ8VHnVyha11dpq8Yh9qxURlCqKSkHxUDFIK3KqqIgBCh44KIKG5N5/zMqQhAw5zcyaWev3uS6u5HmymNyTTPLLOjz3MndHRESkPq3CLkBERHKXQkJERFJSSIiISEoKCRERSUkhISIiKSkkREQkJYWEiIikpJAQEZGUFBIiIpJSm7ALaImOHTt6cXFx2GWIiOSVefPmferu/6cx2+Z1SBQXF1NWVhZ2GSIiecXMPmrstjrcJCIiKSkkREQkJYWEiIikpJAQEZGUFBIiIpKSQkJERFJSSIiISEp5vU5CRLb31PyVjJ+xlFXrN9O1sIDSob0Y0a8o7LIkTTZuqeDe2e9z7oBiOu/eLuOfTyEhEiFPzV/JmCkL2FxRCcDK9ZsZM2UBgIIiz7k7V099l8fmrgCgT9fdGX5g14x/Xh1uEomQ8TOWJgOi2uaKSsbPWBpSRZIOs5asZe8x05MBcfExPbMSEKA9CZFIWbV+c5PmJbet/WILh908Mznu0r4df7tiIAU7tc5aDQoJkQjpWljAynoCoWthQQjVSHNVVTmj/jSPlxavSc49e9lR7N+1fdZr0eEmkQgpHdqLgra1/8osaNua0qG9QqpImmrq/HL2uXp6MiCuPWE/lo87IZSAAO1JiERK9clpXd2Uf1Z89hX/OX5Wcnxgt/Y8+V8DaNs63L/lFRIiETOiX5FCIY9UVFZx6r1v8M+P1yfnZv9yIMUddw2xqm0UEiIiIXnotQ8Z+9dFyfHtp/XllEO6hVjR9hQSIiJZtnj1Ro6/c05yPLh3J+4/p4RWrazB/5vtxZIKCRGRLNlSUcmQO16mfN22K9DmXn0snRq5cjqMxZIKCRGRLLhtxlLunrUsOb7/nBKG9OncpMfY0WJJhYSI5KW495J6a/nnnHbvG8nxGSXdGXfKdzFr+NBSXWEsllRIiEjGxLmX1MYtFRx600t8vbUKgLatjbJrvkf7Xdo2+zHDWCypxXQikjFx7CXl7ox+8h0OvOGFZEBM+skRvHfzsBYFBISzWFJ7EiKSMXHrJfW3JWs4/+Gy5PiSgT258rjeaXv8MBZLKiREJGPi0kuqbiO+ru3bMTNDjfiyvVhSISEiGVM6tFetcxIQrV5S9TXim37Z0fTpunuIVaWXQkJEMibKvaSmzi/n8sf/mRxfe8J+XHj0PiFWlBkKCRHJqKj1ksrVRnyZopAQEWmE+hrxvVw6kL32yI1GfJmikBARaUA+NOLLFIWEiEgKdRvxHdu7E/c1shFfVCgkRETq2FJRybG3v1zr8t2mNOKLEoWEiEgNdRvxPXBuCcfu17RGfFGikBARYftGfGce2p1bfti8RnxRkrGQMLPuwCPAnkAVMMHd7zSzDsDjQDGwHDjd3ddZ4jtxJzAM+AoY6e5vZ6o+ERGADZsrOPTml/gm6LO0U5tWvHXNENoXtKzPUlRkck9iK3CFu79tZrsB88zsRWAkMNPdx5nZaGA0cBVwPLBv8O9w4J7grYjIDjWnHbm7M2bKAv7y1sfJuScuPoJDiztkuty8krGQcPfVwOrg/S/MbDFQBJwEDAw2mwjMJhESJwGPuLsDfzezQjPrEjyOiEi9mtOOvG4jvksH9aR0aPoa8UVJVs5JmFkx0A94E+hc/Yvf3VebWadgsyLg4xr/rTyYU0iISEpNuVtb3UZ8RYUFvPSLYzLSiC8qMh4SZvYt4Eng5+6+cQcnger7gNfzeKOAUQA9evRIV5kikqca0468qsq56JEyZi5Zm5yLWiO+TMloSJhZWxIB8Wd3nxJMr6k+jGRmXYDq71o50L3Gf+8GrKr7mO4+AZgAUFJSsl2IiEi8NNSOvG4jvuuG9+GCo/bOWn35LmMdqYKrlR4AFrv7HTU+NA04N3j/XODpGvPnWEJ/YIPOR4hIQ1Ldre28I4spHv1sMiD6di/kvZuPV0A0USb3JI4EzgYWmNk/grmrgXHAJDO7AFgBnBZ8bDqJy1+XkbgE9rwM1iYiEVG3HXmX9u1o07oVNz27OLlNHBrxZYolLibKTyUlJV5WVtbwhiISCw+8+iE3PrOtEd8dp/flhwfHoxFfU5jZPHcvacy2WnEtInlPjfgyRyEhUkNzFmVJeNSIL/MUEiKB5izKkvCMn7GE3896PzmOeyO+TFFIiASasihLwqNGfNmlkBAJNGZRloRHjfjCoZAQCTS0KEvCoUZ84VJIiARKh/aqdU4CEouySof2CrGqeJu5eA0XTFQjvjApJEQCdRdl6eqm8KgRX+5QSIjUMKJfkUIhRGrEl3sUEiKSE6a8Xc4vJqkRX65RSIhIqD76bBPHjJ+dHPftXsjki4+gbeuM9R+VJlBIiEgoKiqrOPWe1/ln+YbknBrx5R6FhIhknRrx5Q+FhIhkzaJVGxl217ZGfEP268SEs9WIL5cpJEQk47ZUVDL4ttms2rAlOadGfPlBISEiGVU8+tlaYzXiyy8KCRHJiD/MXsatzy9Njofs15n7zjlEjfjyjEJCRNLq48+/4uhbZ9Wam3PlILp32CWkiqQlFBIikhbuzt5jpteau/K4Xlwy8DshVSTpoJAQkRa7ZuoC/vzmilpzy8edEFI1kk4KCRFptndXbmD4716tNfePX32Pwl12CqkiSTeFhIg0WWWV0/Pq2oeWfnNGX07upwVxUaOQEJEm+fH9b/Lqsk+T46LCAl4bPTjEiiSTFBIi0ihz3vuEsx+YW2tuyY3H0a6t7vEQZQoJEdmhLRWV9L7u+Vpz//+Cwzlq344hVSTZpJAQkZQG3DKzViuNo/ftyJ8uODzEiiTbFBIisp26NwACeP/Xw2itRnyxo5AQkaR1m76h340v1pp79rKj2L9r+5AqkrApJEQE2L4R34/79+CmEd8NqRrJFQoJkZj7/axljJ+xtNbch7cMUyM+ARQSEmFPzV/J+BlLWbV+M10LCygd2osR/YrCLitn1NeI79WrBtHt22rEJ9soJCSSnpq/kjFTFrC5ohKAles3M2bKAoDYB0V9jfhGH9+bi4/pGVJFkssUEhJJ42csTQZEtc0VlYyfsTTWITFmygIem6tGfNJ4CgmJpFXrNzdpPurUiE+aSyEhkdS1sICV9QRC18KCEKoJT32N+H57xkGx3puSplFISCSVDu1V65wEQEHb1pQO7RViVdlVtxFf9w4FzLlSjfikaTIWEmb2IDAcWOvuBwRzNwAXAZ8Em13t7tODj40BLgAqgcvcfUamapPoq/5LOY5XN3302SaOGT+71pwa8UlzZXJP4mHgbuCROvO/cffbak6YWR/gTGB/oCvwkpn9h7tXItJMI/oVxSIUqlVUVnHKPa/zTvmG5Jwa8UlLZSwk3P0VMytu5OYnAX9x96+BD81sGXAY8EaGyhOJlPvnfMBNzy5OjnUDIEmXMM5J/NTMzgHKgCvcfR1QBPy9xjblwdx2zGwUMAqgR48eGS5VJLctWrWRYXfNSY6H7NeZCWcfQis14pM0yXZI3APcCHjw9nbgfKC+V7TX9wDuPgGYAFBSUlLvNiJRt/mbSgbfPpvVNdp4z73mWDrt1i7EqiSKshoS7r6m+n0zuw94JhiWA91rbNoNWJXF0kTyxv88v4R7Zr+fHD84soTBvTuHWJFEWVZDwsy6uPvqYHgy8G7w/jTgUTO7g8SJ632BufU8hEhszf3wc07/47bTdGcd1p1fn/xdNeKTjMrkJbCPAQOBjmZWDlwPDDSzg0gcSloO/ATA3Rea2SRgEbAVuFRXNokkbPiqgkNuepGtVYmjqzu3acXca4bQvqBtyJVJHJh7/h7WLykp8bKysrDLEMkId+fKye/wxLzy5Nzki4+gpLhDiFVJFJjZPHcvacy2WnEtkoNeWrSGCx/Z9gfQfw/+Dld8Pz6rxbNBreQbRyEhkkPWbtzCYb+emRwXFRYw84pjtFo6zdRKvvEUEiI5oKrKOX/iW8xe+kly7rmfHc1+XXYPsaroUiv5xlNIiIRs8rxyfvnEP5PjXw3vw/lH7R1iRdGnVvKNp5AQCUndRnwHdS9k8sVH0KZ1q/CKigm1km88hYRIllVUVnHyH17j3ZUbk3OvlA6ixx66t3S2qJV84ykkRLLovlc+4ObpasQXtji3km8qhYRIFixctYET7tp2+9Dv9enMH3+sRnxhilsr+eZSSIhk0OZvKhl42yzWbPw6OadGfJJPFBIiGTLuuSXc+/K2RnwPjTyUQb07hViRSNMpJETS7M0PPuOMCdtuj3LWYT349ckHqBGf5CWFhEiabPiqgoNvepHKoBFfu7aJRny7t1MjPslfCgmRFnJ3Sie/w2Q14pMIUkiItMCLi9ZwkRrxSYQpJESaYc3GLRxeoxFf9w4FvHi5GvFJ9CgkRJqgvkZ8z//8aHrvqUZ8Ek0KCZFGqtuI7/oT+3DekWrEJ9GWMiTMbDpwibsvz145Irln+aebGHjb7OS4X49CnviJGvFJPOxoT+Jh4AUzmwjc6u4V2SlJJDfU14hvzpWD6N6h5Y34dFc0yRcpQ8LdJ5nZs8CvgDIz+xNQVePjd2ShPpFQ3DBtIQ+/vjw5/u0ZB6Xtl7juiib5pKFzEhXAJmBnYDdqhIRIFM1eupaRD72VHH+/T2fuTXMjPt0VTfLJjs5JHAfcAUwDDnb3r7JWlUiWbdxSwYE3vFBr7uXSgey1x65p/1y6K5rkkx3tSVwDnObuC7NVjEgYhv9uTq3zDpm+akl3RZN8sqNzEkdnsxCRbJvydjm/mLTtktY9d2/HG2MGZ7wRn+6KJvlE6yQkduqulgYou3YIHb+1c1Y+v+6KJvlEISGx4e70HfsCG7dsTc797qx+nNi3a9Zr0V3RJF8oJCQW/vjy+9zy3JLk+NDib/PExQNCrCg/aX1H/CgkJNLe/+RLjr395Vpz744dyrd21ku/qbS+I570kyKRVFnl9Lx6eq25P194OEd+p2NIFeU/re+IJ4WERM7Yvy7kodeWJ8cjDurKb8/sF15BEaH1HfGkkJDI+MfH6xnx+9dqzf3rpuPZqY0a8aWD1nfEk0JC8t6Wikp6X/d8rbln/vsoDihqH1JF0aT1HfGkkJC8dvGf5vH8wn8nx/81sCdXHdc7xIqiS+s74kkhIXlp1tK1nFejER/AB78eltZGfLI9re+IH4WE5JX6GvG9UjqIHnu0/B4PIrI9hUQM5euCqLqN+G44sQ8jdftQkYzKWEiY2YPAcGCtux8QzHUAHgeKgeXA6e6+zhId1e4EhgFfASPd/e1M1RZn+bgg6sl55VxR497SXdq3440xx4ZYkUh8ZPLawIeB4+rMjQZmuvu+wMxgDHA8sG/wbxRwTwbrirUdLYjKNWs2bqF49LO1AqLs2iEKCJEsytiehLu/YmbFdaZPAgYG708EZgNXBfOPuLsDfzezQjPr4u6rM1VfXOXDgih358AbXuCLr7c14rv7R/0YfmD2G/GJxF22z0l0rv7F7+6rzaxTMF8EfFxju/JgbruQMLNRJPY26NGjR2arjaBcXxB178vvM65GI77D9u7ApJ8cEWJFIvGWKyeu67tu0evb0N0nABMASkpK6t1GUsvVBVH1NeJbOHYou6oRn0iosv0TuKb6MJKZdQHWBvPlQPca23UDVmW5tljItQVR9TXie/TCwxmgRnwiOSHbITENOBcYF7x9usb8T83sL8DhwAadj8icXFkQdcO0hTz8+vLk+OR+RfzmjIPCK0hEtpPJS2AfI3GSuqOZlQPXkwiHSWZ2AbACOC3YfDqJy1+XkbgE9rxM1SXhm79iHSf/4fVac2rEJ5KbMnl101kpPrTd9YvBVU2XZqoWyQ31NeJ79rKj2L+rGvGJ5CqdFZSsGPVIGS8sWpMcXzKwJ1eqEZ9IzlNISEbNWrKW8x5WIz6RfKWQkIzYsLmCvmPViE8k3ykkJO1OuGsOC1epEZ9IFCgkJG0mzyvnlzX6LHVt347X1WdJJK8pJKTFJr6+nOunLaw1V3btEDp+a+eQKhKRdFFISLO5O72ve56vt1bVmi9o25pX3/s0JxbsiUjLaPWSNMs9s99n7zHTtwsIyN3W4yLSdNqTkCZZtvZLhtzxcoPb5VLrcRFpPu1JSKNUVjnFo5+tFRCPXnQ4RSlajOdK63ERaRmFhDTohmkLa3Vq/WG/IpaPO4EBPTtSOrQXBW1b19o+F1qPi0h66HCTpNSYRny51npcRNJLISHbaWojvlxpPS4i6aeQkFrqNuK7dFBPSoeqEZ9IXCkkBFAjPhGpn0Ii5uprxDfnykF076BGfCKikIi1YXfOYdHqbY34xv5gf84dUBxeQSKScxQSMVS3EV9RYQGvjR4cYkUikqsUEjHy7w1b6H/LzFpzasQnIjuikIgBd+eA62ew6ZvK5Nzvf3QwJxzYJcSqRCQfKCQi7g+zl3Hr89ua7fXfpwN/GXVEiBWJSD5RSERUfY34Fo4dyq4761suIo2n3xgRU1nltfosQaIR34CeHUOqSETymUIiQq5/+l0mvvFRcvzDg4u44/SDQqwovZ6av1I9okSyTCERAW+vWMcPG2jEl++emr+SMVMWsLkicfJ95frNjJmyAEBBIZJBCok81tRGfPls/IylyYCoVn0HPIWESOYoJPLUhRPLeGlxfBrxpbrTne6AJ5JZCok887clazj/4bJac3FoxNe1sICV9QSC7oAnklkKiTwR90Z8pUN71TonAboDnkg2KCTywHG/fYUl//4iOY5jIz7dAU8kHAqJHDap7GOunPxOchz3Rny6A55I9ikkclB9jfjmXTuEPdSIT0SyTCGRIc1Z+OXuXD31XR6buyI5p0Z8IhImhUQGNGfh17K1XzL2rwuZ896ngBrxiUhuUEhkQFMWfn2xpYK7Zr7HQ68tp2Cn1lx/Yh/O7r8XbVpHZ7W0iOSvUELCzJYDXwCVwFZ3LzGzDsDjQDGwHDjd3deFUV9LNWbhl7szdf5KbnluCZ988TVnlHSn9LheugGQiOSUMPckBrn7pzXGo4GZ7j7OzEYH46vCKa1lGlr49e7KDVw/bSHzPlpH327tue+cEg7qXpjtMkVEGpRLxzROAiYG708ERoRYS4uUDu1FQdvWteYK2rbmkoE9uWbqAk68+1WWf7qJW085kKmXHKmAEJGcFdaehAMvmJkDf3T3CUBnd18N4O6rzaxTSLW1WN2FX13at2PAdzoy/oWlfLFlK+ceUczl3/sP2he0DblSEZEdCyskjnT3VUEQvGhmSxr7H81sFDAKoEePHpmqr8WqF36VLf+cXz29kMnzyum/Twdu+MH+9N5z97DLExFplFBCwt1XBW/XmtlU4DBgjZl1CfYiugBrU/zfCcAEgJKSEs9WzU21duMWbnluCVPnr6RL+3bc/aN+nPDdLphFuxGfiERL1kPCzHYFWrn7F8H73wf+HzANOBcYF7x9Otu1pcM3W6t4+PUPufOl96iodC4d1JNLB32HXXbS1cYikn/C+M3VGZga/EXdBnjU3Z83s7eASWZ2AbACOC2E2lrklX99wg1/XcgHn2xicO9O/Gp4H4o77hp2WSIizZb1kHD3D4C+9cx/Bhyb7XrS4ePPv+LGZxbxwqI1FO+xCw+OLGFw785hlyUi0mI6BtICWyoquWf2+9z78vu0MqN0aC8uPHpvdm7TuuH/LCKSBxQSzeDuzFi4hhufWcTK9ZsZfmAXrh62n+6SJiKRo5BoopqN+Hp13o3HLurPET33CLssEZGMUEg0khrxiUgcKSQaULMR36dffs3ph6gRn4jEh0JiB9SIT0TiTiFRj3WbvuG2F5by6NwVdNhlJ2495UBOPaQbrVpptbSIxItCoobKKufRuSu4PWjEN3JAMT8fkhuN+JpzO1QRkZZSSASqG/EtWr2R/vt0YOwPDqDXnruFXRbQvNuhioikQ+xDIh8a8TXldqgiIukU25D4ZmsVD732IXfNzP1GfI25HaqISCbk3m/ELJi/Yh1XPPFPPvhkE8f27sR1Od6Ir6HboYqIZEosV4LtslMbDHhwZAkPjDw0pwMCUt8OtXRor5AqEpG4iOWeRK89d+PFy4/Jm0ta694OVVc3iUi2xDIkgLwJiGrVt0MVEcmmWB5uEhGRxlFIiIhISgoJERFJKXbnJNTeQkSk8WIVEmpvISLSNLE63LSj9hYiIrK9WIWE2luIiDRNrEIiVRsLtbcQEalfrEJC7S1ERJomVieu1d5CRKRpYhUSoPYWIiJNEavDTSIi0jSx25PIZVroJyK5RiGRI7TQT0RykQ435Qgt9BORXKSQyBFa6CciuUghkSO00E9EcpFCIkdooZ+I5CKduM4RWugnIrlIIZFDtNBPRHKNDjeJiEhKORcSZnacmS01s2VmNjrsekRE4iynQsLMWgO/B44H+gBnmVmfcKsSEYmvnAoJ4DBgmbt/4O7fAH8BTgq5JhGR2Mq1kCgCPq4xLg/mREQkBLl2dZPVM+e1NjAbBYwKhl+aWX19KzoCn6a5tnwS9+cP+hqAvgagr0Gq579XYx8g10KiHOheY9wNWFVzA3efAEzY0YOYWZm7l6S/vPwQ9+cP+hqAvgagr0E6nn+uHW56C9jXzPY2s52AM4FpIdckIhJbObUn4e5bzeynwAygNfCguy8MuSwRkdjKqZAAcPfpwPQWPswOD0fFQNyfP+hrAPoagL4GLX7+5u4NbyUiIrGUa+ckREQkh0QqJOLY0sPMupvZLDNbbGYLzexnwXwHM3vRzN4L3n477Fozycxam9l8M3smGO9tZm8Gz//x4EKIyDKzQjObbGZLgtfCETF8DVwe/Ay8a2aPmVm7qL8OzOxBM1trZu/WmKv3+24JdwW/H98xs4Mb8zkiExIxbumxFbjC3fcD+gOXBs97NDDT3fcFZgbjKPsZsLjG+H+A3wTPfx1wQShVZc+dwPPu3hvoS+JrEZvXgJkVAZcBJe5+AIkLX84k+q+Dh4Hj6syl+r4fD+wb/BsF3NOYTxCZkCCmLT3cfbW7vx28/wWJXw5FJJ77xGCzicCIcCrMPDPrBpwA3B+MDRgMTA42ifrz3x34T+ABAHf/xt3XE6PXQKANUGBmbYBdgNVE/HXg7q8An9eZTvV9Pwl4xBP+DhSaWZeGPkeUQiL2LT3MrBjoB7wJdHb31ZAIEqBTeJVl3G+BK4GqYLwHsN7dtwbjqL8W9gE+AR4KDrndb2a7EqPXgLuvBG4DVpAIhw3APOL1OqiW6vverN+RUQqJBlt6RJmZfQt4Evi5u28Mu55sMbPhwFp3n1dzup5No/xaaAMcDNzj7v2ATUT40FJ9guPuJwF7A12BXUkcXqkryq+DhjTr5yJKIdFgS4+oMrO2JALiz+4+JZheU70rGbxdG1Z9GXYk8AMzW07iEONgEnsWhcFhB4j+a6EcKHf3N4PxZBKhEZfXAMAQ4EN3/8TdK4ApwADi9Tqolur73qzfkVEKiVi29AiOvz8ALHb3O2p8aBpwbvD+ucDT2a4tG9x9jLt3c/diEt/zv7n7/wVmAacGm0X2+QO4+7+Bj82sVzB1LLCImLwGAiuA/ma2S/AzUf01iM3roIZU3/dpwDnBVU79gQ3Vh6V2JFKL6cxsGIm/IqtbetwcckkZZ2ZHAXOABWw7Jn81ifMSk4AeJH6ATnP3uie4IsXMBgK/dPfhZrYPiT2LDsB84Mfu/nWY9WWSmR1E4sT9TsAHwHkk/giMzWvAzMYCZ5C44m8+cCGJY+6RfR2Y2WPAQBLdXtcA1wNPUc/3PQjPu0lcDfUVcJ67lzX4OaIUEiIikl5ROtwkIiJpppAQEZGUFBIiIpKSQkJERFJSSIiISEoKCZE0M7OTzczNrHfYtYi0lEJCJP3OAl4lsbhPJK9pnYRIGgU9tJYCg4BpQetukbylPQmR9BpB4r4O/wI+b+yNXURylUJCJL3OItEGguDtWSHWItJiOtwkkiZmtgeJTptrSbRgbh283cv1gyZ5SnsSIulzKok7f+3l7sXu3h34EDgq5LpEmk0hIZI+ZwFT68w9CfwohFpE0kKHm0REJCXtSYiISEoKCRERSUkhISIiKSkkREQkJYWEiIikpJAQEZGUFBIiIpKSQkJERFL6X3DU6a8x6WIZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(df3.A, df3.Y)\n",
    "ax.plot(df3.A, res.predict(df3[['constant', 'A']]))\n",
    "ax.set_xlabel('A')\n",
    "ax.set_ylabel('Y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 11.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Still Program 11.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Let us return to the data in Figure 11.1.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1['constant'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols = sm.OLS(df1['Y'], df1[['constant', 'A']])\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/pettychen/anaconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>   67.5000</td> <td>   19.716</td> <td>    3.424</td> <td> 0.004</td> <td>   25.213</td> <td>  109.787</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>        <td>   78.7500</td> <td>   27.883</td> <td>    2.824</td> <td> 0.014</td> <td>   18.947</td> <td>  138.553</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 11.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 11.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting from the same data as Section 11.2, Program 11.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "df3['A^2'] = df3.A * df3.A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols = sm.OLS(df3.Y, df3[['constant', 'A', 'A^2']])\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/pettychen/anaconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>   -7.4069</td> <td>   31.748</td> <td>   -0.233</td> <td> 0.819</td> <td>  -75.994</td> <td>   61.180</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>        <td>    4.1072</td> <td>    1.531</td> <td>    2.683</td> <td> 0.019</td> <td>    0.800</td> <td>    7.414</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A^2</th>      <td>   -0.0204</td> <td>    0.015</td> <td>   -1.331</td> <td> 0.206</td> <td>   -0.053</td> <td>    0.013</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = res.summary()\n",
    "summary.tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8lFW+x/HPD0IJNZRQEgi9Sjcgiq7oumsXLLiiKArKFnXV3etd9e66u3eruqvrlutrEZGioFhWWMuiYi+UUKSI9JYESCgBQnpy7h8zwYAZCJCZZ+aZ7/v1ymsyT55kflO/ec45zznmnENERKQ6dbwuQEREopdCQkREQlJIiIhISAoJEREJSSEhIiIhKSRERCQkhYSIiISkkBARkZAUEiIiElKC1wWcjtatW7vOnTt7XYaISExZunTpHudcck32jemQ6Ny5MxkZGV6XISISU8xsW033VXOTiIiEpJAQEZGQFBIiIhKSQkJEREJSSIiISEgxPbpJRL7pteVZPDZ/Hdl5haQkJXL/xb0YPTjV67IkRikkRHzkteVZPPjqKgpLywHIyivkwVdXASgo5JSouUnERx6bv+5IQFQqLC3nsfnrPKpIYp1CQsRHsvMKT2q7yIkoJER8JCUp8aS2i5yIQkLER+6/uBeJ9eoetS2xXl3uv7iXRxVJrFPHtYiPVHZOa3ST1BaFhIjPjB6cqlCQWqPmJhERCUlHEiIiMSTSJ0sqJEREYoQXJ0uquUlEJEZ4cbKkjiREJKw0l1Tt8eJkSR1JiEjYVDaPZOUV4vi6eeS15VlelxaTvDhZUiEhImGjuaRqlxcnS6q5SUTCRnNJ1S4vTpZUSIhI2KQkJZJVTSBoLqlTF+mTJdXcJCJho7mkYp+OJEQkbDSXVOxTSIhIWGkuqdim5iYREQlJISEiIiEpJEREJCSFhIiIhKSQEBGRkBQSIiISUthCwsw6mtn7ZrbWzNaY2T3B7S3N7B0z2xC8bBHcbmb2VzPbaGYrzWxIuGoTEZGaCeeRRBnwU+dcH2A4cKeZ9QUeABY453oAC4LXAS4FegS/JgFPhbE2EfGR15ZnMeKP79HlgTcY8cf3NMtsLQpbSDjndjrnlgW/PwSsBVKBUcD04G7TgdHB70cBM1zAQiDJzNqHqz4R8QdNRx5eEemTMLPOwGBgEdDWObcTAkECtAnulgrsqPJrmcFtIiIhxcN05M459uQXszIzj7dW7WTKx5tZsnVfRG477NNymFkT4BXgXufcQTMLuWs121w1f28SgeYo0tLSaqtMEYlRfpiOvKSsgp0HCsnKKyRrfyHZeUVk5wWuV14Wl1Uc9Tt3XtCNoZ1bhr22sIaEmdUjEBDPO+deDW7ebWbtnXM7g81JOcHtmUDHKr/eAcg+9m865yYDkwHS09O/ESIiEl+ifTpy5xwHC8vIzCs46sP/60AoJDe/GHfMp1mbpg1ISUqkT/tmXNS3LSnNG5LaohEpSQ1JTUqkeWK9iNQftpCwwCHDM8Ba59zjVX40DxgP/DF4ObfK9rvM7AXgLOBAZbOUiEgo91/ciwdfXXVUk1MkpyMvK69g96HiIx/4WVWPAILbDpcc3RzWIKEOqUmJpCQlMrJXMqlJX3/4p7ZIpF3zhjRIqBviFiMrnEcSI4CbgVVmtiK47SEC4TDHzCYC24ExwZ+9CVwGbAQKgNvCWJuI+ES4pyPPLy77+sO/ShBk5wWahXYdLKK84ujDgJaN65OS1JCuyY05t0frwId/MBRSWyTSqnF9jtP0HlXMHXuME0PS09NdRkaG12WISIxyzpGbX0zm/qP/88/KKzoSBAcKS4/6nYQ6RvukhqQ0D3zgH/nwr3KZWD86jgJCMbOlzrn0muyr9SREqnhteZYWyPGxgpIyVmYeYPn2PJZv38/yHXnkHio+ap+mDROO/Oef3qkFqS2+/vBPTUokuWkD6taJjaOA2qCQEAmqHG9f2bZdOd4eUFDEoIoKx+Y9h1m+fT8rduSxfHse63YfOtI01KV1Y87r3pr+HZrTqVUjUoJHAs0aRqZDOFYoJESCjjfeXiER/QpKyliydT/LtgWOEFZs38/BojIAmjZIYFBaEnf26cbgtBYM7JhEy8b1Pa44NigkRIL8MN4+njjn2JR7mA/W5fDh+lwWbdlHSVkFdQx6tm3K5QNSGJyWxJC0JLq2bkKdOGoiqk0KCZGgaB9vL3C4uIzPNu09EgyZ+wPPV482TbhleCfO75XM4LQWNGmgj7baokdSJMjr8fbyTc45NuTkHwmFJVv2U1JeQeP6dTmne2t+OLIb5/dMpkOLRl6X6lsKCZGgcI+3l5rJLy7j0417+GBdLh+tzz1ydNerbVNuG9GZ83smk965JfUTtBxOJCgkRKoYPThVoRBhzjnW7T7EB+ty+WBdDhlb91NW4WjSIIER3Vtx14XdOb9nspr9PKKQEJGIO1hUyqcb9vDh+lw+XJ/LzgNFAPRu15Tbz+vKyF7JDElroaOFKKCQEJGwc86xduchPlifwwfrclm2LXC00LRhAuf1aM29PZM5v2cb2jVv6HWpcgyFhIiEhXOOlZkHmLsimzdX7WTXwcDRQt/2zZj0ra6M7NWGwWlJ1Kuro4VoppAQkVq1MSefeSuymPdFNlv3FlC/bh1G9krmJ317MrJnMm2a6WghligkROS0ZecV8vrKbOauyGZN9kHqGJzdrRU/Gtmdi/u1i9jaB1L7FBIickr2Hy7hzdU7mbsim8VbAktpDuyYxMNX9OWKAe11xOATCgkRqbGi0nLmr9nFvBXZfLg+l7IKR7fkxvz0Oz25cmAKnVs39rpEqWUKCRE5oe17C3h+8TbmLNnB/oJS2jdvyMRzu3DVoBT6tm8WMwvoyMlTSIhItcorHO9/lcPMhdv4aEMudcy4qE8bxg3vxIhurTVhXpxQSIjIUXIPFfPiku3MXryDrLxC2jZrwI8v7MHYYWk6jyEOKSREBOcci7bs47mF25i/Zhel5Y4R3Vvx88v7cFHftjqXIY4pJETi2MGiUv61LIvnFm5jQ04+zRomcPPwztw0PI1uyU28Lk+igEJCJA6tyT7Acwu3M3dFFgUl5Qzo0JxHrxvAlQNSSKxf1+vyJIooJETiRFFpOW+u2slzC7exbHseDRLqMGpQCuOGd2JAhySvy5MopZAQ8bltew8za9F25mQEhq92bd2YX1zRl+uGdKB5I50JLcenkBDxofIKx3uVw1fX51K3jvHdvm0ZN7wT53RrpfMapMYUEiI+UlJWwSvLMnnqg01s31dA22YNuPeiHtwwVMNX5dQoJER8oKi0nJcydvDUB5vIPlDEgA7NefDSIRq+KqdNISESwwpLypm1eDuTP9rE7oPFnNmpBb+/pj/n90xWk5LUCoWESAw6XFzGzIXbmPLxZvbklzC8a0ueuH4QZ6u/QWqZQkIkhhwsKmX6p1t55tMt5BWUcl6P1tx9YQ+GdWnpdWniUwoJ8a3Xlmfx2Px1ZOcVkpKUyP0X92L04FSvyzoleQUlTP10K89+uoVDRWVc2LsNd1/YncFpLbwuTXxOISG+9NryLB58dRWFpeUAZOUV8uCrqwBiKij25hcz5ZMtzPx8G/nFZVx8RlvuvrAH/VKbe12axAmFhPjSY/PXHQmISoWl5Tw2f11MhETOwSImf7SZ5xdtp6isnMv6t+fuC7vTu10zr0uTOKOQEF/Kzis8qe3RYueBQv754WZmLd5OWXkFowalcucF3ejepqnXpUmcUkiIL6UkJZJVTSCkJCV6UM2J7dhXwFMfbuLljEwqnOOaIan8aGR3LQcqnlNIiC/df3Gvo/okABLr1eX+i3t5WNU3bd1zmH+8v5F/Lc/CDMakd+SH53ejY8tGXpcmAoQxJMxsKnAFkOOc6xfc9ivgDiA3uNtDzrk3gz97EJgIlAM/ds7ND1dt4n+V/Q7ROrppT34xf357PS8u2U69unUYN7wT3z+/K+2bR+eRjsSvcB5JTAP+Dsw4ZvsTzrk/Vd1gZn2BG4AzgBTgXTPr6ZwrR+QUjR6cGjWhUKm4rJxpn27l7+9tpLC0nFvO7syPLuhGm6aaV0miU9hCwjn3kZl1ruHuo4AXnHPFwBYz2wgMAz4PU3kiEeWc4+0vd/P7N9eybW8BF/Zuw0OX9aF7G63+JtHNiz6Ju8zsFiAD+Klzbj+QCiyssk9mcNs3mNkkYBJAWlpamEsVOX1rdx7kN69/yWeb9tK9TROmTxjG+T2TvS5LpEYiHRJPAb8BXPDyz8AEoLrJZlx1f8A5NxmYDJCenl7tPiLRoGq/Q7PEevzvqDO4cVgaCZqVVWJIREPCObe78nszexp4PXg1E+hYZdcOQHYESxOpNSVlFUz7bAt/WxDodxh/Tmfu/XZPrQInMSmiIWFm7Z1zO4NXrwZWB7+fB8wys8cJdFz3ABZHsjaR0+Wc451gv8PWvQVc0CuZ/7m8r/odJKaFcwjsbGAk0NrMMoFfAiPNbBCBpqStwPcBnHNrzGwO8CVQBtypkU0SS77aFeh3+HRjoN9h2m1DGdmrjddliZw2cy52m/XT09NdRkaG12VIHNubX8zj76xn9uJAv8N9F/XkxrPStBqcRDUzW+qcS6/JvjrjWuQUlJRVMP2zrfx1wYYj/Q73fLsHSY3qe12a1JCfppIPJ4WEyElwzvHu2hx+98aX6neIYX6ZSj4SFBIiNfTVroP89vW1fLJxj/odYlysTyUfSQoJkRM4tt/h11edoX6HGBerU8l7QSEhEkJJWQUzPt/Kkws2UFASmGfp3ovU7+AHsTaVvJcUEiLHcM6xYG0Ov3tzLVv2HA72O/TRwj8+EitTyUcDhYRIFdl5hTw8dzXvrs1Rv4OPRftU8tFEISECVFQ4nlu0jUfe+ooKB/9zWR9uHdFZ/Q4+Fo1TyUcjhYTEvfW7D/HAKytZtj2Pb/VM5nej+2llOJEghYTEreKycv7x3kae+nATTRok8MT3BjJ6UCpm1U1KLBKfFBISl5Zs3ccDr6xkU+5hrh6cys8v70OrJg28Lksk6igkJK4cLCrlkbe+4vlF2+nQIlELAImcgEJC4sZ/Vu/il/NWk3uomNvP7cJPvtuTRvX1FhA5Hr1DxPd2Hyzi4bmrmb9mN33bN+PpW9IZ0CHJ67JEYoJCQnzLOcdLSzP5zb+/pKS8gp9d0pvbz+uiYa0iJ0EhIb6Uc6iIh15dxbtrcxjWpSWPXjuAzq0be12WSMxRSIjvvLFyJz9/bRWHS8r5+eV9mDCiC3XqaFiryKlQSIhv5BWU8PDcNcz7IpuBHZrz5+sHar4lkdOkkBBfeP+rHH72ykr2HS7hp9/pyQ9HdiMhivsetCqaxAqFhMS0/OIyfvfGl8xevINebZsy9dah9Ett7nVZx6VV0SSWhPxXy8zeNLPOkStF5OQs3LyXS/7yES8u2cEPzu/GvLtHRH1AwPFXRROJNsc7kpgGvG1m04FHnXOlkSlJ5PiKgh+oUz/dQlrLRrz0g7M5s1NLr8uqMa2KJrEkZEg45+aY2RvAw0CGmc0EKqr8/PEI1CdylC925PGTOSvYlHuYm4d34sHLesfcWdNaFU1iyYl69kqBw0ADoOkxXyIRU1Hh+L8PNnLNU59RUFLOzInD+M3ofjEXEBBYFS2xXt2jtmlVNIlWId9hZnYJ8DgwDxjinCuIWFUiVezNL+Ync77gw/W5XD6gPb+/uj/NE+t5XdYp06poEkuO92/Y/wBjnHNrIlWMyLEWb9nH3bOXsb+glN+O7sdNZ6X5Yr0HrYomseJ4fRLnRbIQkaoqKhxPfbiJP7+9jk6tGjP11qGckRL9I5f8Tud3xJ/Ya9AV39uTX8x9L67g4w17uHJgCr+/uh9NG8Zu85Jf6PyO+KSQkKiycPNefjx7OXmFpfz+6v6MHdbRF81LfnC88zsUEv6lkJCoUDl66fF31tO5VWOm3TaMvinNvC5LqtD5HfFJISGeq9q8dNXAFH5/TX+aNNBLM9ro/I74FL0zoElc+HzTXi578mMWb9nHH67pz5M3DFJARCmd3xGf9G4UT5RXOP7x/kb+8m6geWn6hGH0aa/mpWim8zvik0JCIi73UKB56ZONexg9KIXfXq3mpVih8zvij96ZccjLse6fbdrDPS+s4GBhKY9c25/r0zV6SSSaha1PwsymmlmOma2usq2lmb1jZhuCly2C283M/mpmG81spZkNCVdd8a5yrHtWXiGOr8e6v7Y8K6y3W17hePLdDYybsohmDROYe9cIvjfUH2dPi/hZODuupwGXHLPtAWCBc64HsCB4HeBSoEfwaxLwVBjrimterGWQc6iIm59ZxBPvrmfUoFTm3XUuvdup/0EkFoStuck591E1ixaNAkYGv58OfAD8LLh9hnPOAQvNLMnM2jvndoarvngV6bHun23cw49fWEF+cSmPXjuAMekddPQgEkMiPQS2beUHf/CyTXB7KrCjyn6ZwW3fYGaTzCzDzDJyc3PDWqwfhRrTXttj3csrHE+8s56bnllE88QE5t55LtcPVf+DSKyJlvMkqvvkcNXt6Jyb7JxLd86lJycnh7ks/4nEWPecQ0WMm7KIJxds4Opg81KvdlqCRCQWRXp00+7KZiQzaw/kBLdnAh2r7NcByI5wbXEh3GPdP9mwh3tfXE5+cRmPXjeAMWeqeUkklkU6JOYB44E/Bi/nVtl+l5m9AJwFHFB/RPiEY6x7eYXjyQUb+Nt7G+iW3IRZdwynZ1sdPYjEurCFhJnNJtBJ3drMMoFfEgiHOWY2EdgOjAnu/iZwGbARKABuC1ddUvv25hdz16zlfL55L9cO6cBvRp8Rk8uKisg3hXN009gQP/p2Nfs64M5w1SLhs3bnQW6fnsGe/GIeu24AY9I7nviXRCRm6N89OWXz1+zivhdX0LRhAnO+fzYDOyZ5XZKI1DKFhJw05wKT8/3p7fUM7NCcybek07ZZQ6/LEpEwUEjISSkqLef+l1fy7y+yGTUohUeuHUDDY4bUioh/KCSkxnYdKGLSzAxWZR3gvy/pxQ/P76bhrSI+p5CQGlmxI49JMzI4XFzG5JvT+U7ftkd+5uWssiISXgoJOaG5K7K4/+WVtGnagBkTzzlqcr7KWWUrJw2snFUWUFCI+EC0TMshUaiiwvHIf77inhdWMKhjEnPvHPGN2Vu9mFVWRCJHRxJSrfziMu59YQXvrt3N2GEd+fVV/aif8M3/KSI9q6yIRJZCQr5hx74Cbp+ewcbcfH51ZV/Gn9M5ZAd1SlIiWdUEQm3PKisi3lBzkxxl4ea9XPX3T9h5oJBptw3l1hFdjjuCKRKzyoqId3QkIUfMWrSdh+euJq1VI6bckk7X5CYn/J1wzyorIt5SSAhl5RX89o21TPtsK9/qmczfxg6meWK9Gv9+OGaVFZHooJCIcwcKSrlz1jI+2biHCSO68NBlvUmoq1ZIEQlQSMSxjTn53DEjg8z9BTx67QCuH6oZXEXkaAqJOPXxhlx+9Pwy6tetw6w7hjO0c0uvSxKRKKSQiEOvLM3kZ6+spHubJkwZn06HFo28LklEopRCIo4453jqw008+p91nN21Ff+85UyaNax5B7WIxB+FRJwor3D8+t9rmPH5Nq4cmMKfxgygQYKm+BaR41NIxIGi0nLue3EFb63exR3ndeHBS/tQp46m+BaRE1NI+NyBglLumJHB4q37+Pnlfbj9vK5elyQiMUQh4WNZeYXcOnUx2/YW8Lexg7lyYIrXJYlIjFFI+NRXuw5y69QlHC4uY9qEoZzTrbXXJYlIDFJI+NDnm/YyaUYGjRrU5aUfnv2NNSBilVbAE4k8hYTP/PuLbH465ws6tWrEtAnDSPXJlN1aAU/EG5qkx0ee+WQLd89ezsCOzXnpB2f7JiBAK+CJeEVHEj5QUeH4w1trefrjLVxyRjv+csMgGtbz1zkQWgFPxBsKiRhXXFbO/S+tZN4X2dxydid+eeUZ1PXhORBaAU/EG2puimGHikq57dklzPsim59d0ptfX+XPgACtgCfiFR1JxKjdB4u49dklbNh9iMevH8g1Qzp4XVJYaQU8EW8oJGLQxpx8xk9dTF5BCVNvHcq3eiZ7XVJEaAU8kchTSMSYpdv2MXF6Bgl16vDi98+mX2pzr0sSER9TSMSQt9fs4u7Zy0lJSmT6bcNIa6V1IEQkvBQSYVLbZwfPydjBA6+sZECHJKbeOpSWjevXYrUiItVTSIRBbZ8dPPWTLfzv619yXo/W/PPmM2lUX0+biESGhsCGQW2dHeyc44l31vO/r3/Jpf3aMWV8ugJCRCLKk08cM9sKHALKgTLnXLqZtQReBDoDW4HrnXP7vajvdNXG2cEVFY7fvPElz366lTFnduAP1/Qnoa4yXUQiy8tPnQucc4Occ+nB6w8AC5xzPYAFwesxKdRZwDU9O7isvIL/fmUlz366lQkjuvDItQMUECLiiWj65BkFTA9+Px0Y7WEtp+V0zg4uLivnrlnLeXlpJvdd1JNfXKGlRkXEO141cDvgbTNzwD+dc5OBts65nQDOuZ1m1qa6XzSzScAkgLS0tEjVe1JO9ezgotJyJs1cykfrc3n4ir5MOLdLJMoVEQnJnHORv1GzFOdcdjAI3gHuBuY555Kq7LPfOdfieH8nPT3dZWRkhLnayCgoKWPitAwWbtnLH6/pz/eGRmcAikjsM7OlVZr6j8uT5ibnXHbwMgf4FzAM2G1m7QGClzle1OaFQ0WljJ+6mEVb9vL49QMVECISNSIeEmbW2MyaVn4PfBdYDcwDxgd3Gw/MjXRtXjhQUMq4ZxazfHsefxs7hKsH+3uiPhGJLV70SbQF/mVmlbc/yzn3HzNbAswxs4nAdmCMB7VF1P7DJYx7ZhHrdx/i/24awnfPaOd1SSIiR4l4SDjnNgMDq9m+F/h2pOvxSu6hYm5+ZhFb9hxm8i3pXNCr2n56ERFP6fRdD+w+WMSNTy8kO6+IqbcOZUT31l6XJCJSLYVEhGXlFXLj0wvZc6iY6ROGMaxLS69LEhEJSSERQTv2FTD26YUcKCxl5u1nMSTtuCN8RUQ8p5CIkM25+dw0ZRGFpeXMun04/TtosSARiX4KiQjYsPsQN05ZREWFY/Ydw+nTvpnXJYmI1IhCIsy+zD7IuGcWkVDHePH7w+nepqnXJYmI1JhCIoxWZx3gpimLaFS/LrPuGE6X1o1P+W/V9kp3IiI1oZAIk692BY4gmjRI4IVJw+nY8tTXo67tle5ERGoqmqYK942NOfmMm7KIhgl1mX3H6QUE1N5KdyIiJ0shUcu27T3MTVMWAsbzd5xFWqvTCwionZXuREROhUKiFgVOlFtESVkFz99+Ft2Sm9TK3z3dle5ERE6VQqKWVE61caiolJkTz6JXu9obxXQ6K92JiJwOdVzXgj35xdz49EL25pcwc+Iw+qXW7olyp7rSnYjI6VJInKb9h0sYN2UR2XlFTJ8wjMFhmmpj9OBUhYKIRJxC4jQcKCzllqmL2bznMM/eOlST9YmI76hP4hTlF5dx27OL+WrXQf457kxN9y0ivqQjiVNQWFLOxGlL+CLzAP+4cQgX9NaCQSLiT3EXEqc7vUVRaTmTZmaweOs+/vK9QVzST0uOioh/xVVInO70FiVlFdw1axkfb9jDY9cNYNQgdSSLiL/FVZ/E6UxvUVZewT0vLOfdtTn8dnQ/xqR3DFeZIiJRI65C4lSntyivcPzXS1/w1upd/OKKvowb3ikc5YmIRJ24ColTmd6iosLx0KureG1FNvdf3IuJ53YJV3kiIlEnrkLiZKe3cM7xq3+v4cWMHfz4wu7ceUH3SJQpIhI14qrj+mSmt3DO8Ye3vmLG59uY9K2u3PednpEuV0TEc3EVElDz6S2eeGc9kz/azPizO/Hgpb0xswhUJyISXeIuJGriH+9v5K/vbeSGoR355ZVnRCwgtESpiEQbhcQxpny8mcfmr2P0oBR+d3V/6tSJXEBoiVIRiTZx1XF9IjMXbuO3b6zlsv7t+NOYgdSNUECAligVkeikkAiak7GDX7y2mov6tOEv3xtMQt3IPjRaolREopFCApi7IoufvbKS83q05u83DqF+QuQfFi1RKiLRKO5D4j+rd/KTOV9wVpeWTL45nYbHnEcRKVqiVESiUVx3XL/31W7unr2cgR2a88z4oSTW9yYgQEuUikh0ituQ+GTDHn7w3DJ6t2vGtAnDaNzA+4dCS5SKSLSJy+amxVv2cfuMJXRt3ZgZE4bRrGE9r0sSEYlKURcSZnaJma0zs41m9kA4biO5aQOGdm7Jc7efRYvG9cNxEyIivuB9G0sVZlYX+AfwHSATWGJm85xzX9bm7XRp3ZiZE8+qzT8pIuJL0XYkMQzY6Jzb7JwrAV4ARnlck4hI3Iq2kEgFdlS5nhncJiIiHoi2kKhuHgx31A5mk8wsw8wycnNzI1SWiEh8iraQyASqLh7dAciuuoNzbrJzLt05l56cnBzR4kRE4k20hcQSoIeZdTGz+sANwDyPaxIRiVtRNbrJOVdmZncB84G6wFTn3BqPyxIRiVtRFRIAzrk3gTe9rkNERKKvuUlERKKIOedOvFeUMrNcYFs1P2oN7IlwOdEk3u8/6DEAPQagxyDU/e/knKvRyJ+YDolQzCzDOZfudR1eiff7D3oMQI8B6DGojfuv5iYREQlJISEiIiH5NSQme12Ax+L9/oMeA9BjAHoMTvv++7JPQkREaodfjyRERKQW+CokIrFgUbQxs45m9r6ZrTWzNWZ2T3B7SzN7x8w2BC9beF1rOJlZXTNbbmavB693MbNFwfv/YnCaF98ysyQze9nMvgq+Fs6Ow9fAfcH3wGozm21mDf3+OjCzqWaWY2arq2yr9nm3gL8GPx9XmtmQmtyGb0KiyoJFlwJ9gbFm1tfbqiKiDPipc64PMBy4M3i/HwAWOOd6AAuC1/3sHmBtleuPAE8E7/9+YKInVUXOk8B/nHO9gYEEHou4eQ2YWSrwYyDdOdePwLQ+N+D/18E04JJjtoV63i8FegS/JgFP1eQGfBMSxOmCRc65nc65ZcHvDxH4cEglcN+nB3ebDoz2psLwM7MOwOXAlOB1Ay4EXg7u4vf73wz4FvAMgHOuxDmXRxy9BoISgES4XPMuAAADCklEQVQzSwAaATvx+evAOfcRsO+YzaGe91HADBewEEgys/Ynug0/hUTcL1hkZp2BwcAioK1zbicEggRo411lYfcX4L+BiuD1VkCec64seN3vr4WuQC7wbLDJbYqZNSaOXgPOuSzgT8B2AuFwAFhKfL0OKoV63k/pM9JPIXHCBYv8zMyaAK8A9zrnDnpdT6SY2RVAjnNuadXN1ezq59dCAjAEeMo5Nxg4jI+blqoTbHcfBXQBUoDGBJpXjuXn18GJnNL7wk8hccIFi/zKzOoRCIjnnXOvBjfvrjyUDF7meFVfmI0ArjKzrQSaGC8kcGSRFGx2AP+/FjKBTOfcouD1lwmERry8BgAuArY453Kdc6XAq8A5xNfroFKo5/2UPiP9FBJxuWBRsP39GWCtc+7xKj+aB4wPfj8emBvp2iLBOfegc66Dc64zgef8PefcTcD7wHXB3Xx7/wGcc7uAHWbWK7jp28CXxMlrIGg7MNzMGgXfE5WPQdy8DqoI9bzPA24JjnIaDhyobJY6Hl+dTGdmlxH4L7JywaLfeVxS2JnZucDHwCq+bpN/iEC/xBwgjcAbaIxz7tgOLl8xs5HAfznnrjCzrgSOLFoCy4FxzrliL+sLJzMbRKDjvj6wGbiNwD+BcfMaMLNfA98jMOJvOXA7gTZ3374OzGw2MJLAbK+7gV8Cr1HN8x4Mz78TGA1VANzmnMs44W34KSRERKR2+am5SUREaplCQkREQlJIiIhISAoJEREJSSEhIiIhKSREapmZXW1mzsx6e12LyOlSSIjUvrHAJwRO7hOJaTpPQqQWBefQWgdcAMwLTt0tErN0JCFSu0YTWNdhPbCvpgu7iEQrhYRI7RpLYBoIgpdjPaxF5LSpuUmklphZKwIzbeYQmIK5bvCyk9MbTWKUjiREas91BFb+6uSc6+yc6whsAc71uC6RU6aQEKk9Y4F/HbPtFeBGD2oRqRVqbhIRkZB0JCEiIiEpJEREJCSFhIiIhKSQEBGRkBQSIiISkkJCRERCUkiIiEhICgkREQnp/wH4CPj2kJZu/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df3.sort_values('A', inplace=True)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(df3.A, df3.Y)\n",
    "ax.plot(df3.A, res.predict(df3[['constant', 'A', 'A^2']]))\n",
    "ax.set_xlabel('A')\n",
    "ax.set_ylabel('Y');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = df3.shape[0]\n",
    "yvar = (res.resid * res.resid).sum() / (n - 2)  # = res.mse_resid\n",
    "xval = np.array([[1, 90, 90 * 90]])\n",
    "X = df3[['constant', 'A', 'A^2']]\n",
    "XpXinv = np.linalg.inv(np.dot(X.T, X))\n",
    "se_mean = np.sqrt(yvar * np.dot(xval, np.dot(XpXinv, xval.T)))[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24.2462945968487"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "se_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate      95% C.I.\n",
      "E[Y|A=90]   197.13   (144.75, 249.51)\n"
     ]
    }
   ],
   "source": [
    "t = scipy.stats.t.ppf(0.975, n - 3)\n",
    "ypred = res.predict(xval)[0]\n",
    "print('           estimate      95% C.I.')\n",
    "print(\n",
    "    'E[Y|A=90]   {:>6.2f}   ({:>6.2f}, {:>6.2f})'.format(\n",
    "        ypred, ypred - t * se_mean, ypred + t * se_mean\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Note: This confidence interval is slightly different from what we get in the Stata code."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}