P2CH3_02 Linear Regression Assumptions on Residuals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Residuals related assumptions\n",
    "\n",
    "In this notebook, we work on the linear regression assumptions related to the residuals of a model.\n",
    "\n",
    "The residuals of a reliable linear regression model should be\n",
    "\n",
    "- normally distributed\n",
    "- not auto-correlated\n",
    "- and with a constant variance\n",
    "\n",
    "We work on the Sales ~ sqrtTV + Radio + TV:Newspaper model, which is the best we've designed so far.\n",
    "\n",
    "Let's import the packages and load the advertising dataset into a pandas dataframe\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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>TV</th>\n",
       "      <th>Radio</th>\n",
       "      <th>Newspaper</th>\n",
       "      <th>Sales</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>230.1</td>\n",
       "      <td>37.8</td>\n",
       "      <td>69.2</td>\n",
       "      <td>22.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>44.5</td>\n",
       "      <td>39.3</td>\n",
       "      <td>45.1</td>\n",
       "      <td>10.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>17.2</td>\n",
       "      <td>45.9</td>\n",
       "      <td>69.3</td>\n",
       "      <td>9.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>151.5</td>\n",
       "      <td>41.3</td>\n",
       "      <td>58.5</td>\n",
       "      <td>18.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>180.8</td>\n",
       "      <td>10.8</td>\n",
       "      <td>58.4</td>\n",
       "      <td>12.9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      TV  Radio  Newspaper  Sales\n",
       "0  230.1   37.8       69.2   22.1\n",
       "1   44.5   39.3       45.1   10.4\n",
       "2   17.2   45.9       69.3    9.3\n",
       "3  151.5   41.3       58.5   18.5\n",
       "4  180.8   10.8       58.4   12.9"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.stats.api as sms\n",
    "\n",
    "\n",
    "df = pd.read_csv('./data/advertising.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run the regression and get the following result summary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>Sales</td>      <th>  R-squared:         </th> <td>   0.931</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.930</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   876.7</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Thu, 18 Jul 2019</td> <th>  Prob (F-statistic):</th> <td>2.89e-113</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>16:03:14</td>     <th>  Log-Likelihood:    </th> <td> -346.84</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   200</td>      <th>  AIC:               </th> <td>   701.7</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   196</td>      <th>  BIC:               </th> <td>   714.9</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>          <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>    <td>   -1.3430</td> <td>    0.346</td> <td>   -3.880</td> <td> 0.000</td> <td>   -2.026</td> <td>   -0.660</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>np.sqrt(TV)</th>  <td>    0.9368</td> <td>    0.029</td> <td>   31.832</td> <td> 0.000</td> <td>    0.879</td> <td>    0.995</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Radio</th>        <td>    0.1902</td> <td>    0.007</td> <td>   27.553</td> <td> 0.000</td> <td>    0.177</td> <td>    0.204</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>TV:Newspaper</th> <td> 5.644e-05</td> <td> 2.58e-05</td> <td>    2.187</td> <td> 0.030</td> <td> 5.55e-06</td> <td>    0.000</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 5.590</td> <th>  Durbin-Watson:     </th> <td>   2.064</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.061</td> <th>  Jarque-Bera (JB):  </th> <td>   6.074</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.250</td> <th>  Prob(JB):          </th> <td>  0.0480</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.692</td> <th>  Cond. No.          </th> <td>2.37e+04</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 2.37e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                  Sales   R-squared:                       0.931\n",
       "Model:                            OLS   Adj. R-squared:                  0.930\n",
       "Method:                 Least Squares   F-statistic:                     876.7\n",
       "Date:                Thu, 18 Jul 2019   Prob (F-statistic):          2.89e-113\n",
       "Time:                        16:03:14   Log-Likelihood:                -346.84\n",
       "No. Observations:                 200   AIC:                             701.7\n",
       "Df Residuals:                     196   BIC:                             714.9\n",
       "Df Model:                           3                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "================================================================================\n",
       "                   coef    std err          t      P>|t|      [0.025      0.975]\n",
       "--------------------------------------------------------------------------------\n",
       "Intercept       -1.3430      0.346     -3.880      0.000      -2.026      -0.660\n",
       "np.sqrt(TV)      0.9368      0.029     31.832      0.000       0.879       0.995\n",
       "Radio            0.1902      0.007     27.553      0.000       0.177       0.204\n",
       "TV:Newspaper  5.644e-05   2.58e-05      2.187      0.030    5.55e-06       0.000\n",
       "==============================================================================\n",
       "Omnibus:                        5.590   Durbin-Watson:                   2.064\n",
       "Prob(Omnibus):                  0.061   Jarque-Bera (JB):                6.074\n",
       "Skew:                          -0.250   Prob(JB):                       0.0480\n",
       "Kurtosis:                       3.692   Cond. No.                     2.37e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 2.37e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = smf.ols('Sales ~ np.sqrt(TV) + Radio + TV:Newspaper' , data = df).fit()\n",
    "result.summary()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normality of the residuals\n",
    "\n",
    "\n",
    "Let's first check the normality of the residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1166104a8>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXycV33v8c9vRjPaV2u1VsvyJjteZTuJ44QkZCMBh7JkgwQKmLSkQCkt3EJ59V5a2suF0hYCbpqmNIGQBEKaFBycAEnATrzItux4iR1J1m5b+y7NaGbO/WPGQVEkeySN9Mzye79eelkzz3NGv7Glr4/Oc55zxBiDUkqp6GWzugCllFJzS4NeKaWinAa9UkpFOQ16pZSKchr0SikV5eKsLmAy2dnZpqyszOoylFIqYhw8eLDTGJMz2bGwDPqysjKqq6utLkMppSKGiDROdUyHbpRSKspp0CulVJTToFdKqSinQa+UUlFOg14ppaKcBr1SSkU5DXqllIpyGvRKKRXlNOiVUirKheWdsUqp6Xl8X9Os2t+9uSRElahwpD16pZSKchr0SikV5TTolVIqygUV9CJys4icEpFaEfnyJMe3ichREakRkWoRuSrYtkoppebWJS/GiogdeBC4AWgBDojIc8aYE+NO+w3wnDHGiMhq4ClgeZBtlVLM/oKqUlMJpke/Cag1xtQbY9zAE8C28ScYYwaNMSbwMBkwwbZVSik1t4KZXlkINI973AJsnniSiLwf+AcgF7h1Om0D7bcD2wFKSnSql1LjdQ+52Xemi7qOQRx2G4kOOysXprG2OBO7TawuT4W5YHr0k30XmXc8YcwzxpjlwO3A16fTNtD+IWNMlTGmKidn0t2wlIo5Y14fP61u5tsvnGJPbScJcXbsInQMuHj6UCv/+ps3OXWu3+oyVZgLpkffAhSPe1wEtE11sjHmdyKyWESyp9tWKfUHI24vj+1toKFrmK0V2VxZkU16ogMAYwzH2/p58cR5Hn2tkQ9uKGJdSabFFatwFUzQHwCWiMgioBW4E7h7/AkiUgHUBS7GrgecQBfQe6m2Sql3GnJ5eHh3PZ2Dbu7cWMzqooy3HRcRVhWmszQvlcf2NvCzgy0AGvZqUpcMemOMR0QeAHYBduARY8xxEbk/cHwH8AHgXhEZA0aAOwIXZydtO0fvRamoYIzh6UMtdA66+diVZSzOSZnyXGecjY9eXsajgbBPTXBQkTv1+So2yR8my4SPqqoqU11dbXUZSs2rC9Mr99R28svXz3Lb6gKuXJwdVFu3x8eDL9Xi9vr43PVLSHDYp/W1da2byCciB40xVZMd0ztjlQojrb0j/OrYOVbkp3JF+YKg2znjbHxwQxH9I2PsfP3sHFaoIpEGvVJhwmcMzxxuITnezgfWFyEyvWmTxVlJbF2SQ3VjD6fPD8xRlSoSadArFSaONPfS1jvKzavySYqf2Qri16/IJSclnl8cPYsvDIdllTU06JUKAyNuLy+cOE9hRuI7ZthMh8Nu492VeXQOujja0hfCClUk06BXKgz8x+56+kbGeM9lBdimOWQz0cqFaeSnJfDbN9q1V68ADXqlLNc77OYHL9dRWZDGouzkWb+eTYTrlucGevW9IahQRToNeqUs9qO9jQy5vbx7RV7IXrNSe/VqHA16pSw0Oublh6828K5lOeSnJ4TsdW0iXLMsh85BN7XtgyF7XRWZNOiVstCFO2A/ffXikL/2yoVpJDvt7D/THfLXVpFFg14pi3h9hn//XT1ritK5vDwr5K8fZ7OxoTSTN8710z86FvLXV5FDg14pi7x44hwNXcN8+prF0745KlhVZVn4DBxq7JmT11eRQYNeKYs8treRwoxEblqZP2dfIzslnvLsZA40dOtF2RimQa+UBeo7BtlT28Vdm4rnfIeojYuy6Bkeo04vysYsDXqlLPCT/U3E2YQPVxVf+uRZWlmQRoLDxhGdUx+zNOiVmmejY15+drCFG1fmkZsWuimVU4mz26gsSOPE2X48Pt+cfz0VfjTolZpnvzp2jp7hMe7ZXDpvX3PVwnRGx3zUtQ/N29dU4UODXql59vi+JhZlJ09rvfnZqshNIT7OxrE2XegsFmnQKzWPmrqG2d/QzQc3FGGb44uw48XZbawoSONEWz9en86+iTUa9ErNo2cOtyICt68rnPevvWphOiNjXuo7dfZNrNGgV2qemMAOUpcvWkBhRuK8f/0leSk442wca9Xhm1ijQa/UPDnc3EtD1zDvXz//vXnwb0qyPD+VE2cH9OapGDOz/cqUUpN6fF/TlMeerWnFYRcGRz0XPW8uLctL5WhLH2d7RynMnP/fKpQ1tEev1DzweH0cbeljRUEaCQ67ZXUsyUsF4HS7bh4eSzTolZoHp84PMDLmZV1xpqV1pMTHUZiRyOlzGvSxJKigF5GbReSUiNSKyJcnOX6PiBwNfLwqImvGHWsQkddFpEZEqkNZvFKR4nBTLynxcVTkplhdCkvzUmjqHmbE7bW6FDVPLhn0ImIHHgRuASqBu0SkcsJpZ4BrjDGrga8DD004fq0xZq0xpioENSsVUYbdHk6dG2BNUfqcL2AWjKV5qRigtkOnWcaKYHr0m4BaY0y9McYNPAFsG3+CMeZVY8yFBa/3AkWhLVOpyHW0pQ+vMawrsXbY5oKizCQSHDYdvokhwQR9IdA87nFL4LmpfAJ4ftxjA7wgIgdFZPv0S1QqstU095KbGk9BCPeEnQ27TViSm8rp9gGMTrOMCcEE/WS/a0763SEi1+IP+i+Ne3qLMWY9/qGfz4jI1VO03S4i1SJS3dHREURZSoW/rkEXTd3DrC/JnLNdpGZiaV4KA6MezvWPWl2KmgfBBH0LMH7R7CKgbeJJIrIaeBjYZozpuvC8MaYt8Gc78Az+oaB3MMY8ZIypMsZU5eTkBP8OlApjh5t7EWBNcYbVpbxNRa5/mmVdh65mGQuCCfoDwBIRWSQiTuBO4LnxJ4hICfBz4KPGmNPjnk8WkdQLnwM3AsdCVbxS4cwYQ01zL+U5yaQnOqwu523SEx0sSHZyRi/IxoRL3hlrjPGIyAPALsAOPGKMOS4i9weO7wC+BiwAvh/49dQTmGGTBzwTeC4OeNwY86s5eSdKhZmm7mG6h9xctzzX6lImtSg7mWNtfbocQgwIagkEY8xOYOeE53aM+/yTwCcnaVcPrJn4vFKx4FBTLw67sHJhmtWlTKo8J4Xqxh7O9uk4fbTTO2OVmgNjXh+vt/aycmE68XHWLXlwMeXZyYB/o3IV3TTolZoDp84NMDrmY22YXYQdLy3RQXaKk3q9IBv1NOiVmgOHm3pITYhjcY71Sx5czKLsFBq6hvB4ddPwaKZBr1SIDbo8nDo/wJqijLBY8uBiynOScXl8HG/rt7oUNYc06JUKsaMtvfgMrA+TJQ8uZlFgnH5vfdclzlSRTINeqRA71NjDwowE8sNkyYOLSUtwkJ0Sz74z3VaXouaQBr1SIXS2b4S2vtGI6M1fULogiUNNPfh8Op8+WmnQKxVChxp7sIuwpih8Z9tMVJqVRO/wGPWdOvsmWmnQKxUiY14fNc29LC9IJTk+crZjLlmQBPj/k1LRSYNeqRB55VQHQ25vRA3bAGSnxJOR5OCgBn3U0qBXKkR+drCF5Pg4lgY24I4UNhHWl2RS3agXZKOVBr1SIdA95OY3b5xnXXH4z52fzIbSTOo6hugZcltdipoDGvRKhcBzNa2MeQ3rSiLnIux4G0r9w02Hm3X4Jhpp0CsVAj871MLKhWkUpCdaXcqMXLiLV8fpo5MGvVKz9Ma5fo619vPBDUVWlzJjiU47KxemadBHKQ16pWbp6YMtOOzCtrWFVpcyKxtKM6lp7mVMFziLOhr0Ss2C2+Pj54dauW55LlnJTqvLmZUNpZmMjvk4eVYXOIs2kXNXh1Lz4PF9TdM6/2hLL11DbgrSE6fdNtxcuCB7sLGH1RF0Z6+6NO3RKzUL+890k5nkoCI3vNedD0ZBeiIL0xN0nD4KadArNUOdAy7qO4fYWJaFTSJv7vxkNpRl6VIIUUiDXqkZ2t/QjU3+MOQRDTaUZNDWN0pb74jVpagQ0qBXagY8Xh+HmnpYUZBGaoLD6nJCZkNpFoAO30QZDXqlZuD11j6G3V42LcqyupSQWl6QSqLDrkEfZTTolZqB1+q7yE6JpyLMN/+eLofdxtriDA41adBHEw16paapuXuYlp4RrijPQqLkIux4G0ozOd7Wz7DbY3UpKkSCCnoRuVlETolIrYh8eZLj94jI0cDHqyKyJti2SkWavfVdOONsrIuwdeeDtaE0E6/PcLSlz+pSVIhcMuhFxA48CNwCVAJ3iUjlhNPOANcYY1YDXwcemkZbpSLGoMvD0dY+1pdkkuCwW13OnLiwAqeO00ePYHr0m4BaY0y9McYNPAFsG3+CMeZVY8yF74q9QFGwbZWKJPvPdOP1GS4vj66LsONlJDlZnJOs8+mjSDBBXwg0j3vcEnhuKp8Anp9uWxHZLiLVIlLd0dERRFlKzS+P18e++i6W5KaQm5pgdTlzan1JJoebezHGWF2KCoFggn6yq02T/uuLyLX4g/5L021rjHnIGFNljKnKyckJoiyl5teRll4GXB6uWpJtdSlzbn1pJt1Dbhq7hq0uRYVAMEHfAhSPe1wEtE08SURWAw8D24wxXdNpq1S4M8bw+zc7yU9LiLoplZO5sMG5TrOMDsEE/QFgiYgsEhEncCfw3PgTRKQE+DnwUWPM6em0VSoSvNk+SPuAi6uWZEfllMqJKnJTSImP06CPEpdcptgY4xGRB4BdgB14xBhzXETuDxzfAXwNWAB8P/BD4AkMw0zado7ei1JzZvebnaQlxLG6KN3qUuaF3Sb+G6cae60uRYVAUOvRG2N2AjsnPLdj3OefBD4ZbFulIklb7wi1HYPctDKfOFvs3GO4viSD771Uy5DLQ3K8bl0RyWLnu1apGdpT24nTbmNTWfROqZzMutJMfMZ/EVpFNg16pS6ib2SMIy29VJVlkuiMzhukprK+2H9B9nCTBn2k06BX6iJeq+vEGLhycfRPqZwoPcnB4pxkDusF2YinQa/UFFxjXvY3dLOyMD3iN/6eqfUlmRxq0hunIp0GvVJTqG7sYXTMx9aK2OvNX7CuRG+cigYa9EpNwusz7KnrpHRBEsVZSVaXY5n1pf4FznQ+fWTToFdqEsfb+ugdHmNrRWwvx7EkN1VvnIoCGvRKTWCMYXdtJwuSnSwvSLW6HEvpjVPRQYNeqQkauvw7SG2pyMYWA8sdXMr6kgzeONfPkEt3nIpUGvRKTbC7tpMkp/2thb1i3YUbp3THqcilQa/UOJ0DLt4428/mRQtwxumPB8C6Yr0gG+n0O1mpcXbXdWK3SVTvIDVdGUlOyvXGqYimQa9UQNegi0ONPawtziA1wWF1OWFFb5yKbBr0SgX8aG8THp/hqhi+QWoq6/XGqYimQa8UMDrm5bG9DSzLSyU3Lbr3g52JDaX+C9PVumF4RNKgVwp47kgbnYPumNgPdiaW5KaQnuiguqHb6lLUDGjQq5hnjOGHe/y9+fLsZKvLCUs2m1BVmsl+DfqIpEGvYt7+M92cONvPx7eUxcR+sDO1cVEW9R1DdA26rC5FTZMGvYp5P3y1gYwkB9vWFlpdSljbWOYfpz/QoOP0kUaDXsW01t4Rdh0/x50bS2JuB6npWlWYTnycTcfpI5AGvYppj73WCMBHryi1uJLwFx9nZ01xBgc06COOBr2KWSNuLz/Z38RNK/MpzEi0upyIsLEsk2Nt/Qy7dYGzSKJBr2LWf9e00jcyxse3LLK6lIixsSwLr8/ohuERRoNexaQLUyorC9LeusioLm19aSYi/plKKnIEFfQicrOInBKRWhH58iTHl4vIayLiEpEvTjjWICKvi0iNiFSHqnClZuO1+i5OnR/gYzqlclrSEhxUFqSx70yX1aWoabhk0IuIHXgQuAWoBO4SkcoJp3UDnwW+NcXLXGuMWWuMqZpNsUqFyg/3NJCV7OR9axZaXUrEuaJ8AYeaehkd81pdigpSMD36TUCtMabeGOMGngC2jT/BGNNujDkAjM1BjUqF1Nm+EX598jx3biwmwaFTKqfrisULcHt8uj59BAkm6AuB5nGPWwLPBcsAL4jIQRHZPtVJIrJdRKpFpLqjo2MaL6/U9PysugWfgTs3llhdSkTauCgLm8DeOh2+iRTBBP1kA5jTWZR6izFmPf6hn8+IyNWTnWSMecgYU2WMqcrJyZnGyysVPJ/P8GR1M1sqFlCyIMnqciJSWoKDywrTea1egz5SxAVxTgtQPO5xEdAW7BcwxrQF/mwXkWfwDwX9bjpFKhUqr9Z10dIzwl/etMzqUsLK4/uapnV+eqKDPbVd/HBPAx/bUjY3RamQCaZHfwBYIiKLRMQJ3Ak8F8yLi0iyiKRe+By4ETg202KVmq0nq5tJT3Rw08p8q0uJaOU5KXiNobFryOpSVBAu2aM3xnhE5AFgF2AHHjHGHBeR+wPHd4hIPlANpAE+Efk8/hk62cAzgelrccDjxphfzc1bUerieobc7Dp2jrs3l+hF2FkqXZCETaC+U4M+EgQzdIMxZiewc8JzO8Z9fg7/kM5E/cCa2RSoVKg8c7gVt9fHHRuLL32yuqj4ODtFmUnUdwxaXYoKgt4Zq2KCMYYnDzSzpiidFQVpVpcTFRbnJNPaO0L/qM6qDnca9ComHGnp49T5Ae7QKZUhU5Gbis/Aq7WdVpeiLiGooRulIslkM0ieOdyCwy6MjnmnPcNETa4kK4n4OBuvnO7g5lUFVpejLkJ79CrquTxejrT0cVlhhl6EDSG7TajITeGVUx0YM51ba9R806BXUe9Yax9uj09XqZwDS3NTaesbpbZdL8qGMw16FfUONPSQkxJPSZbeCRtqS/JSAHjltC5bEs406FVUa+8fpal7mKqyTF2OeA5kJDn9wzca9GFNg15FterGHmwC60p02GauXLM0h3313bq9YBjToFdRy+PzL6W7oiCNlHidYDZXrlmag9vrY68ucha2NOhV1Dp5doBht5eNZVlWlxLVNi3KIslp5zcn260uRU1Bg15FrYON3aQnOqjITbG6lKiW4LBzzdIcXjxxHp9Pp1mGIw16FZV6h928eX6QDaWZ2PQi7Jy7cWUe7QMualp6rS5FTUKDXkWlg43+be426EXYeXHdsjzibMILx89bXYqahAa9ijo+YzjY2MPi3BQyk51WlxMT0pMcXF6+gBdOnLO6FDUJDXoVderaB+kdGaOqVHvz8+mmlXnUdwxR2z5gdSlqAg16FXUONPaQ5LRTqcsRz6sbKv27du3S4Zuwo0Gvokr3kJuTbf2sK84gzq7f3vMpPz2BNcUZ7DquwzfhRn8SVFT5+aEWvMZQpXPnLXHrZfkcbenjjG4xGFY06FXUuLCLVHFmInlpCVaXE5Peu2YhIvBcTZvVpahxNOhV1DjU1Mub7YN6J6yFCtIT2bwoi2drWnWN+jCiQa+ixlMHmkly2rmsMN3qUmLa7WsLqe8c4lhrv9WlqAANehUVBl0e/udoG+9dvZB43UXKUresKsBpt/HfNa1Wl6ICNOhVVPjFkTaG3V7u2FRsdSkxLz3JwbuW5fA/R9rw6to3YUGDXkWFJ6ubWZqXwrriDKtLUcC2tYW0D7jYU9tpdSmKIINeRG4WkVMiUisiX57k+HIReU1EXCLyxem0VWq2Tp0b4HBTLx+uKtZdpMLEuytzyUxy8MSBJqtLUcAld2MQETvwIHAD0AIcEJHnjDEnxp3WDXwWuH0GbZWalScPNOOwC3+0vsjqUmLS4/smD/OVC9P51bFz/NsrdaQmOCY95+7NJXNZmgoIpke/Cag1xtQbY9zAE8C28ScYY9qNMQeAsem2VWo2XB4vzxxu4caV+WTpAmZhZWNZFj4DhwIriSrrBBP0hUDzuMctgeeCEXRbEdkuItUiUt3RoRsNq+C8eOI8PcNj3FGlF2HDTU5qPIuykznQ2INP59RbKpign2zQM9h/taDbGmMeMsZUGWOqcnJygnx5FeuePNBMYUYiV1VkW12KmsTGsiy6h9zUd+iSCFYKJuhbgPHdpSIg2PubZ9NWqYtq7h5md20nH64qxmbTi7DhaOXCNJKcdt043GLBBP0BYImILBIRJ3An8FyQrz+btkpd1JMHmhHgg1V6ETZcOew2NpVlcfJsP12DLqvLiVmXDHpjjAd4ANgFnASeMsYcF5H7ReR+ABHJF5EW4AvAV0WkRUTSpmo7V29GxQ63x8cTB5q4bnkuhRmJVpejLuLyxQuw2YQ9dTqn3iqXnF4JYIzZCeyc8NyOcZ+fwz8sE1RbpWZr1/FzdA66uefyUqtLUZeQluBgTVEGBxt7ePeKPJKcQcWOCiG9M1ZFpB/tbaQ4K5FrluiF+0hwVUU2Y17D/jPdVpcSkzToVcR58/wA+850c/emUr0IGyHy0xNYkpvCa3VdjHl9VpcTczToVcT58b4mnHYbH9aLsBHl6qU5DLg8HGjQXv1806BXEWXY7eHpgy3cclk+C1LirS5HTcPinBQWZSfzyqkO7dXPMw16FVGeq2ljwOXhI3oRNiJdvyKXAZdHx+rnmQa9ihjGGB7b28jy/FSqSjOtLkfNQHl2CuXZybxyugO3R3v180WDXkWMmuZejrf1c8/lpboccQS7fkUegy4Pr+ndsvNGg15FjB/tbSLZaef964JdU0+Fo0XZySzPT+XlU+106t2y80KDXkWEniE3vzjaxu3rCkmJ1xtuIt0tqwoY8/r4pxdPW11KTNCfGBWWJm5m8dKpdlweHwtS4qfc6EJFjpzUeDaXL+CJ/U3ce0Upy/PTrC4pqmmPXoU9j9fH3rouluSmkJ+WYHU5KkSuX55LaoKD//M/JzC6Xv2c0qBXYe9ISx8DLo+uOR9lkpxxfPGmZbxa18XTh1qtLieqadCrsGaMYU9tJ3lp8VTkplhdjgqxezaVUFWayd/98oRemJ1DGvQqrNV2DHKuf5SrKrJ1SmUUstmEf/ijyxhyefj6L05YXU7U0qBXYe2V0x2kxsexpijD6lLUHFmSl8qfvquCZ2va2HX8nNXlRCUNehW2mrqHqe8Y4qol2cTZ9Vs1mn3m2gpWFabxpaePcq5v1Opyoo7+9Kiw9fKpdhIddjYtyrK6FDXHnHE2/uXOdbjGfPz5kzV4fToLJ5Q06FVYausd4Y1zA2ypWEB8nN3qctQ8WJyTwv9+30peq+/iBy/XWl1OVNGgV2Hp5dMdxMfZuKJcp1TGkg9VFfG+NQv59oun+e0b560uJ2po0Kuwc/JsP8db+7iifAGJTu3NxxIR4f9+YDWVBWl87ic11LYPWl1SVNCgV2HnW7tOEe+wsVX3g41JiU47D91bRbzDxqceraZL59fPmq51o8LKwcZufvNGOzdW5mlvPoYVZiTybx/dwN3/vo+P//AAj3/q8ncsZjfbNY/u3lwyq/aRRHv0KmwYY/jmr06RneLkysU6Nh/rNpRm8f171nO8rZ/tj1bj8nitLiliadCrsPHy6Q72nenmgWsrcMbpt6byb1LyzQ+s5tW6Lu5/7CCjYxr2MxHUT5OI3Cwip0SkVkS+PMlxEZF/DRw/KiLrxx1rEJHXRaRGRKpDWbyKHmNeH3/3ixOULUjirhj6lVpd2gc2FPGN91/GS6c6+NSj1Rr2M3DJoBcRO/AgcAtQCdwlIpUTTrsFWBL42A78YMLxa40xa40xVbMvWUWjx15rpK5jiK/eWqnz5tU73L25hG9+cDW7azu575H99I2MWV1SRAnmYuwmoNYYUw8gIk8A24DxKxBtAx41/kWl94pIhogUGGPOhrxiFXW6h9z8869Ps3VJNtevyLW6HDWPpntB9UMbinn6YAs3fucV7ruijIwk5xxVFl2CGbopBJrHPW4JPBfsOQZ4QUQOisj2mRaqote3XjjFkNvL126r1BUq1UWtLc7gvivL6B0eY8crdbouTpCCCfrJfvImLkRxsXO2GGPW4x/e+YyIXD3pFxHZLiLVIlLd0dERRFkqGlQ3dPP4viY+dmUZS/JSrS5HRYCK3BS2X10OwL/9ro66Dr2p6lKCCfoWoHjc4yKgLdhzjDEX/mwHnsE/FPQOxpiHjDFVxpiqnBy9USYWuDxevvT0UQozEvnCDUutLkdFkIL0RO6/ZjHpiQ5+uKeBI829VpcU1oIJ+gPAEhFZJCJO4E7guQnnPAfcG5h9cznQZ4w5KyLJIpIKICLJwI3AsRDWryLYgy/VUdcxxN+/fxXJ8XrvnpqejCQnn756McVZSTxZ3czv3+zQvWencMmfLmOMR0QeAHYBduARY8xxEbk/cHwHsBN4D1ALDAMfDzTPA54JjLvGAY8bY34V8nehIs7xtj5+8HItt69dyLuW6QVYNTOJTjsf31LGTw+28Pyxc/SNjPGeywqw6bWetwmqG2WM2Yk/zMc/t2Pc5wb4zCTt6oE1s6xRRZkRt5fPPVFDZpKTr713pdXlqAjnsNu4c2MxzyfEsaeui/6RMT5UVYxDN6t5i/5NqHn3jZ0nqW0f5NsfXkNWsk6PU7NnE+HW1Qt5z6p8jrX18597zjDi1hurLtCgV/PqxRPneWxvI5/aukhXp1Qhd9WSHO7cWExzzwiP7DnDsNtjdUlhQYNezZsznUN84akaVhWm8cWbllldjopSq4sy+MjmEs73j/Ifu88w7NKw16BX82LI5WH7o9XE2YQdH9mgyxyoObUsP42PXF5Kx4CLh3efYTDGw16DXs05n8/wxZ8eoa5jkO/etZ6izCSrS1IxYGleKvdeUUbnoIuHf1/PwGjsro+jk5fVnBi/hsnO18+yu7aTW1bl09Q9POsNI5QKVkVuCvddWcajrzXw8O4zfGpr+Ts2MIkF2qNXc2pPbSe7azu5vHwBV1XoZiJq/i3OSQmsj+OO2dk4GvRqztQ097Dz9bNUFqRx2+oCXbBMWaY8O4V7NpfS3u/iv15rwO3xWV3SvNKgV3OiprmXn1a3sCg7mTs2FuudispyS/NSuWNjMc3dw/xoX2NMbU2oQa9C7tmaVn5a3UxZdjL3XlGmdyiqsLGqMJ0/Wl9Ibfsgn/3JYXq2ZDgAAAoLSURBVDze2OjZ60+gCqn/erWBzz9ZQ1l2MvddUaZ7v6qws6E0i1svK2DX8fP81dNH8fmifyG02Lv8rOaEMYZ/evE03/1tLTdU5nFVRbb25FXY2lKRzdK8VL7z69Okxsfxt+9bGdXXkDTo1awNuz385U+P8svXz3LnxmL+7vZVPFXdYnVZSl3UZ6+vYGB0jId3nyE1wRHVd2tr0KtZae4eZvtjBzl1rp+/fs9yPrW1PKp7Rip6iAhfuXUFgy4P33upltSEOD59zWKry5oTGvRqxp5//Sx/9fRRAP7jYxu5VteVVxFGRPj791/GgMvDPzz/BqkJDu7eXGJ1WSGnQa+mbcjl4Rs7T/LjfU2sKUrnu3etp2SBLmugIpPdJnznw2sZcnn4yn+/TnK8nW1rC60uK6Q06NWUJluqoL5zkKcPttA7PMbWimxuWJnH7tpO/95iSkUoZ5yNHR/ZwH2P7OfPn6zB7fHxoariSzeMEBr0KijDLg+7TpzjQEMPWclOPrm1nEXZyVaXpVTIJDjs/OfHN/Lpxw7ylz87ysiYl3uvKLO6rJDQoFcX5TOGgw097DpxjtExL1dVZPPuFXk6P15FpSRnHP9+bxUPPH6Yrz17nPP9o/zFDcuw2SJ7goEGvZqUMYbT5wd4/thZzve7KF2QxLY1heSnJ1hdmlJzKsFh5wcfWc/Xnj3Ogy/VcaZziG9/aC2JzsjdQ0GDXr2Nz2f49cnzPPhyHUeae8lKdnL3phJWLkzTaZMqZjjsNr7x/lUszknm73eepL5jD9+9ax1L8lKtLm1GNOgVAB6vj1++fpbvv1THqfMDlGQlcfvaQtaXZhBn02EaFXtEhE9uLaciN4W/eOoI7/3ebv7mtkru2lgScUM5GvTzYLYbbczlvN72gVF+Wt3CT/Y30dIzwpLcFP75jrXctrpA725VCnjXslye/9xWvvDUEb7yzDF+fqiVr29bReXCNKtLC5oGfQwaGB3j1yfP82xNG7vf7MTjM1xRvoCv3lrJjZV5EddbUWqu5aYl8Ogfb+LpQy384/NvcNt3f8/tawv502sXU5Eb/sM5GvSzNOTycLZvhNbeUdp6R2jrHaFjwIXL48Pt8eHy+GjqHgLAGWcnIc5GfJyNeIedZKed9EQH6YlO0pMcJDvtczIOPub1caKtn/1nunnpVDsHGroZ8xoKMxL5xNZF3FFVTHlOSsi/rlLRxGYTPlRVzI2V+Xz3t2/y431NPFPTynXLcvnghiKuW5EbtpveBxX0InIz8C+AHXjYGPOPE45L4Ph7gGHgY8aYQ8G0DWden6FjwEVrIMAvfLwV6n0j9A6/fcNhu01YkOwkwWHHYReccXaGA1uX9QyN4fJ4cQX+A5jIbpNA8Ps/0hLiSE3wf56fHk96opPkeDtJjjgSnfa3pjgaYxgd89E74qZ7yE1rzwgtPSPUdgzyxtl+Tp4dYGTMX8OyvFQ+cVU5N1Tmsq44U3vvSk1TepKDr95WyZ9eW8F/7jnDkwea+c0b7aQlxLGlIputS3KoKsukPDuZuDBZwVWMufhazCJiB04DNwAtwAHgLmPMiXHnvAf4M/xBvxn4F2PM5mDaTqaqqspUV1dP+814vD48PoPPGLy+cR/G4PPB6JiXQZeHIZeHYfcfPh90eegfGaNj0EXHgIuOQTedAy7O94/imbBWdVpCHAszEinMSGThWx8Jbz3OTY1/xz/uZGP0PmMYcnnoGxmjf2SM3pEx+sZ99I+M0T/qwXuRtbLjbIKBKc9JT3SwPD+VFQVpbCzLoqosk7y04KdH6ibeKpqF6tqX12fYU9vJL4+e5XdvdnC2bxSABIeNitwUFqb7syE/PYGC9AQWJMeT6LST5LST7Ax02uw24uxCnF1m/FuBiBw0xlRNdiyYHv0moNYYUx94sSeAbcD4sN4GPGr8/2vsFZEMESkAyoJoGzKr/nYXo2Mz2zFGBBYkx5OT6v+oyEkhLy2ewszEt4K9ID2B1ARHSGq1iZCa4PC/Xubk5xhjGHZ72VKRzfn+UfpHxxh2exlyeRhxexke8yL4p4IlOOxkJDnITHK8VW9WslOnRCo1x+w24eqlOVy9NAdjDHUdQxxt6eVYaz+1HYM0dA3xWl0XAy7PJV8rO8VJ9VdvCHmNwQR9IdA87nEL/l77pc4pDLItACKyHdgeeDgoIqeCqC2kGqY+lA10zlcd80TfU/iLtvcDYfSe7gndS4XsPTUC8jczbl461YFggn6yLuHEsYKpzgmmrf9JYx4CHgqinnknItVT/UoUqfQ9hb9oez+g78kqwQR9CzB+GbcioC3Ic5xBtFVKKTWHgrkkfABYIiKLRMQJ3Ak8N+Gc54B7xe9yoM8YczbItkoppebQJXv0xhiPiDwA7MI/RfIRY8xxEbk/cHwHsBP/jJta/NMrP36xtnPyTuZWWA4pzZK+p/AXbe8H9D1Z4pLTK5VSSkW28JjNr5RSas5o0CulVJTToA+SiPyZiJwSkeMi8k2r6wkVEfmiiBgRyba6ltkSkf8nIm+IyFEReUZEMqyuaaZE5ObA91utiHzZ6npmS0SKReQlETkZ+Bn6nNU1hYqI2EXksIj8wupapqJBHwQRuRb/Hb2rjTErgW9ZXFJIiEgx/uUpomWtgxeBVcaY1fiX3vhfFtczI4GlQx4EbgEqgbtEpNLaqmbNA/yFMWYFcDnwmSh4Txd8DjhpdREXo0EfnD8B/tEY4wIwxrRbXE+ofAf4K6a4iS3SGGNeMMZcuM98L/77NiLRW8uOGGPcwIWlQyKWMebshYUOjTED+IOx0NqqZk9EioBbgYetruViNOiDsxTYKiL7ROQVEdlodUGzJSLvA1qNMUesrmWO/DHwvNVFzNBUS4pEBREpA9YB+6ytJCT+GX9naWaLbM0TXY8+QER+DeRPcugr+P+eMvH/yrkReEpEyk2Yz029xHv6a+DG+a1o9i72nowxzwbO+Qr+oYIfz2dtIRT00iGRRkRSgKeBzxtj+q2uZzZE5Dag3RhzUETeZXU9F6NBH2CMefdUx0TkT4CfB4J9v4j48C9k1DFf9c3EVO9JRC4DFgFHAqtbFgGHRGSTMebcPJY4bRf7dwIQkfuA24Drw/0/4osIZtmRiCMiDvwh/2NjzM+tricEtgDvCyzTngCkiciPjDEfsbiud9AbpoIQuAt4oTHmayKyFPgNUBLBQfI2ItIAVBljwmJVwZkKbHLzT8A1xpiw/k/4YkQkDv/F5OuBVvxLidwdoXeVA29tTvRfQLcx5vNW1xNqgR79F40xt1ldy2R0jD44jwDlInIM/4Wx+6Il5KPM94BU4EURqRGRHVYXNBOBC8oXlg45CTwVySEfsAX4KHBd4N+mJtATVvNAe/RKKRXltEevlFJRToNeKaWinAa9UkpFOQ16pZSKchr0SikV5TTolVIqymnQK6VUlPv/u059svO+9zoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(result.resid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Jarque Bera statistic tests\n",
    "that the residuals have the skewness and kurtosis of a Normal distribution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The JB test for our model has a pvalue of 0.048, very close to the 0.05 threshold.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.04797124511944102"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sms.jarque_bera(result.resid)[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we cannot reject the null hypothesis that the distribution of the residuals\n",
    "is shaped like a normal distribution.\n",
    "\n",
    "To confirm that let's run a Kolmogorov Smirnoff test on the residuals\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.08135122894405838, pvalue=0.13407824555999018)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "stats.kstest(result.resid , 'norm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-value of 0.134 also means that we **cannot** reject the null hypothesis\n",
    "that the residuals are normally distributed.\n",
    "\n",
    "We reach the same conclusion with both tests\n",
    "and it's fair to assume that the residuals are close to being normally distributed\n",
    "\n",
    "## Auto correlation of the residuals\n",
    "\n",
    "The next assumption we would like to check is that\n",
    "the residuals are not correlated.\n",
    "\n",
    "We can do that with the Durbin-Watson test\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.063722393242939"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from statsmodels.stats.stattools import durbin_watson\n",
    "durbin_watson(result.resid, axis=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which gives a value of 2.06\n",
    "\n",
    "indicating that the residuals are not auto correlated\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## heteroskedasticity\n",
    "\n",
    "Finally we can check  the heteroskedasticity assumption\n",
    "This assumption means that the\n",
    "spread of the residuals should stay constant.\n",
    "\n",
    "If we plot the residuals vs fitted values plot for our model we get\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD5CAYAAAA+0W6bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dbZBU13kn8P8zQ8P0WJRa2sE2NEjIKhckRJYoTRGnoGSDVuBVHGkiraNis462vBvyIVZFrHY2yKmyUZwtUTUba11JKtmJnIpd1rLIQh5LkbJIG0goUZG0oBmQKcESEcBuWIkuMWPwjETPzNkP0z303Dnn3nNv39e+/1+VSky/3Hv6wjz39HOec44opUBERNnTkXQDiIgoGAZwIqKMYgAnIsooBnAiooxiACciyigGcCKijFrg9QIRWQHgewA+CWAawKBS6tsishPAbwO4WH/p15RSL7sdq6enR61cubKlBmfZ5OQkFizwvOS5w+uix+uil8frcuTIkapSaonzcZurMAngMaXUWyKyGMAREXm1/txTSqn/atuIlStX4vDhw7YvbzvVahU9PT1JNyN1eF30eF308nhdROSs7nHPAK6UugDgQv3Pl0XkHQDlcJtHRER++foeIiIrAawF8AaA9QC+KiK/BeAwZnrplzTv2QZgGwAsX74c1Wq1xSZn19jYWNJNSCVeFz1eFz1el2usA7iIXAdgL4BHlVI/E5E/B/BNAKr+/z8G8BXn+5RSgwAGAaC3t1fl7auPU94/vwmvix6vix6vywyrKhQRKWAmeD+jlHoeAJRS7ymlppRS0wD+EsC66JpJREROngFcRATAdwC8o5T6VtPjS5te9usAfhx+84iIyMQmhbIewJcBvC0iI/XHvgZgq4jcgZkUyhkAvxNJC4mIEjA0XMHAvpM4PzqBZaUi+resQt/adNVv2FShvAZANE+51nwTEWXV0HAFjz//NiZqUwCAyugEHn/+bQBIVRDnTEwiIoeBfSdng3fDRG0KA/tOJtQivXxNZyKiXGg1/XF+dMLX40lhD5yI2koj/VEZnYDCtfTH0HDF+hjLSkVfjyeFAZyI2koY6Y/+LatQLHTOeaxY6ET/llWhtDEsTKEQUVsJI/3RSLdkvgqFiCgpzlz2xtVLcOrsBbx5oWYMqstKRVQ0wdpv+qNvbTl1AduJKRSiNjI0XMH6Xftxy46XsH7Xfl9537TR5bK///o5VK9cdc1tZyX9EQb2wInaRFZql23pctlOzbnt5p76g3eWceDExVSnP8LAAE7UJtwG77IYvGxz1o0bVfONa++RCp584LZMfm4/mEIhahNZqV22ZZuz7hTJxKSbKDCAE7WJrNQu29Llsp2KhU5MKaV9Lqs3Lj+YQiFqE/1bVs1JJQDJD97pZkQCduV5ulK+RhXK6cvXqlAG9p0MpeokixjAidpE2mqXdYOq/c8dBRRQm1azj7kNtOpK+arVpfM2dEjbjSsuDOBEbSRNtcu6QdXa1Px0R6sDrWm7ccWJAZyIIuEnB91qvjpNN644MYATUSRMMyJNr82aNGz4wCoUIoqEroqk0CkodMzdHyaL+eowVjwMAwM4EUWib20ZTz5wG8qlIgRAuVTEwL++HQNfuh3leo+7UcM9sO9kpqb9p2XDB6ZQiCgybrlpXYXKzheOY2zCvFBVXLzSI2mZNMUeOBHFzlShMjpRSzQlAdilR9IyaYoBnIhiZ9NTTWo6vE16JC0rHjKAE1HsbHuqSUyHt0mP6PL7SSyexRw4UU4lWQanm/avk0R5oe2GEGmoPWcPnCiHbMvgotogwtmDvaG7kJrywrSkR2ywB06UQzZrh0e9QYSzB6vbPm1g30ls3zMS6zeELE3NZwAnyiGbPK/XYF7YAa45oLvdPDasWNTSefy2Jc2YQiHKIZsyOFOQbwTTKGchpmWiTNoxgBPlkC7PKwA2rl4y+7MpyHcIIg+uaZkok3ZMoRDlgK7i5ME7y3jm9XNoLPCqAOw9UkHvzTeib21ZWylS6BTtkrBA8OA6NFzBEy8ex6XxGgCgVCzg+mIBoxO1ea/N4qJXUWIAJ2pDzQG71F3AlQ8n522i0FXogDMUNw9kNnLAO184PhtMJ6f1wRsIFlyHhivof+7onJvC6EQNHQIUOmS2zUB6K0GSxABOlDFe9dvOAcBGz7bZRG3KWIPt7El/NDk9+2fD9pMAECi4Duw7qe3RTyvg+uICdC9cMO9zVqtV3+dpVwzgRBliU9qnGwD0o7knbXusUrEQqGrDLe0yOl7D8Nc3+z5mnnAQkyhDbKozbHPRpWLBc8KKzbGKhU7svG+N1Tmd3NIuzHd7YwAnyhCb6gybwNcIul7reXgdSwA8eGfwmun+LatQ6JR5jxc6hPluC0yhEGWIzTodpuqRjy1coF1r2y34eq1ZogAcOHEx4Ke5dm5nFcrO+9ZkYiJN0jwDuIisAPA9AJ8EMA1gUCn1bRG5EcAeACsBnAHwG0qpS9E1lYh0AdWZ9vCaCt4YBLWZot58LNP+ltyQODk2PfBJAI8ppd4SkcUAjojIqwD+HYC/U0rtEpEdAHYA+P3omkqUfa2uAGi7TocpKAZZ36RxrPW79lut0kfx8QzgSqkLAC7U/3xZRN4BUAZwP4DP11/2XQB/DwZwIqOwFodqpcdqs4hVo63Om4RN75/i5SsHLiIrAawF8AaAT9SDO5RSF0Tk44b3bAOwDQCWL1+e6xrOsbGxpJuQSnm5LrsPHseyrkmgq/nRSew+eFy7QFMU16Vrahy3LtY8MTU++7v52qkqBg+eRtfUND61eOa5P3l5GPes+TievHcldr95DtUrV9Fz3UJsXXcTNqxYFOvvdV7+vdiwDuAich2AvQAeVUr9TGT+yLGOUmoQwCAA9Pb2qp6eniDtbBt5//wmebgub16oQWH+7827l2u4/+lj2lRImNdlaLiCM1c6MKWZjVMuFWfPNfD0MVRGFeBo6+nXL+Kph8rY83vJ12bn4d+LDasyQhEpYCZ4P6OUer7+8HsisrT+/FIA70fTRKL24JYrdlvRL4xNFRrpG13wBoDxq5OzxzUNSiqAqwGmjGcAl5mu9ncAvKOU+lbTUy8AeLj+54cB/Cj85hG1D90KgM10K/rZ7pzjxWtG5aXxGvqfO4qh4YrrjYarAaaLTQ98PYAvA9gkIiP1/+4FsAvAPSJyCsA99Z+JyKB5GzGTyujEbE/7tVPV0NbFtgm8tSmFJ148jv4tqzSJnhmsOEkXmyqU1+BMhl1zd7jNIWpvXiV5AGZ72oMHT9dz0fP57QmbJgA5XRqvoW9tGYfPfjBnqVmAFSdpxKn0RAnwSqcAM6sAdhqKBfz2hG3O1+yP+m7DUw/d4TrNnpLHqfRECXBOyDGt0jqlFIqFznm11xtXL8H6XfutJwQ5zwdAe85SsTDnPQzY6SbKbYHfkPX29qrDhw/Hdr60qVarLH/S4HWBNqVy62KF05cFv/nZm3DgxEXXFEix0Omrhzw0XEH/D47O2TChoVy/IQDp3Jk9j/9eROSIUqrX+ThTKEQpYBo4bCwW1bxXpY7fgc2+tWUMfOn22QHV5nNXRifQ/9xR9P/gaKQbF1PrGMCJUqBvbdmYRjk/OoHdb/zE8xh+Bzb71pZxaMcmlEvFeeeuTal5vXPuCp8+zIETRcDvolVDwxUI9Hlp2wqSZaVioMWy/AR+1oGnC3vgRCELMvlmYN9JbfAWzKRXTNUoDY2BzSCTfvxUtLAOPF0YwIlCFmTyjdv09b61ZWz95RXG9zZK/A6cuBho0o+uxLDQKSh0zL1psA48fZhCIQpZkI0PTGmSxiDjH/Xdhn++eAWH3v1g9rn1t96IZ377V2Z/3r5nxPd5AfMa47rH0lCFQtcwgBOFyCuXbaJba3vRgg48snnV7HHfOjd3GdW3zo1haLgyG1RttlszMdV8M2CnGwM45UYru+HYvtcrl22i6wVvu2vZnMe9NmLYuHoJvv/6uXnH9ipBpOxiAKdcaGU3HNN7D5/9AAdOXJwT1N1y2dv3jGBg30lj8Hf2gps3STAdtzI6gfW79qN/yyrj5sKtbDpM6cZBTMqFVlb1M733mdfPzav4KHUX9AcBWpoQ47WW+PY9I5FtOkzpxQBOmWez4YEpiNkEN7dedbOJ2hRGx2uexwsyIcZrMSq3BTFY+te+GMAp02xrrk1BzCa4+QmAtisLBZk16bWWODB/3WeW/rU3BnDKNNvUiK4HaxvcdO+12xEWoS0HC8yd+m6iAC4BmyMcxKRMc0uNOCtHHryzPG/Q0Sa46SpENq5egr1HKq7blAHm5WD7t6xyrWwZGq5g98HjePNCbd5z/VtWYfueEW1vv1wq4tCOTZ6fidoDAzhlmqn2+fpiYV7lyN4jlcA9Ul2ddO/NN84G4A4R427v/VtWzQv+T7x4HJea8uXNVTEA8Pjzb2NZ1yQUZF7FDHfMoQauBx6jPK5jbKOV6+Is8QNmAllXoWNOgGyIqodqaofzhqF7nbN9wExAv3WxwruXZc5zzW1v9OAroxPorN9AGv8vt/HMyTz+HnE9cGpLzYN7zXlfUzVIVCV1pnY4A6jX7vDnRyesK2b61pZn8/ON3n/j/1y/Ox+YQqHM06U3Gj1TJz+Dh35nbtpsQeZ1A1nW1AM3PdfM7YbgnKlJ7YcBnFKplWnvgH5tET85Yt3sy+17RvDonpGWthxzW9u70b5Gflv3nJPXDYGTeNobAziljs20d68Ab1phz/YmoOvZNkaLKqMT6P/BUUBmdq4xtVFHd2MBZjYT3nnfGgDA3iOVOYOTAuDBO/W9e6/NHjiJp70xgFPqeNV273zhOEYn9BUcziAeNH3g1XPVbQZsk7LwurGs37Vfe+M4cOKi9qZluiEArErJAwZwSh23hZtM9c9h5HubA6QIEKRAy9l20zcFUzvdPrvuW8mTD9yGJx+4bV41SjtXodA1DOCUOm5pAbeY2kq+15m2CVpd25yyCLICoumzd4oYv5Uc2rGJgTqnWEZIqRN06nor+V6v8j4bzpRFkBUQTVP+dZOEAA5S5h0DOKWOrqbaq0PcPD3da2VCnSCBsFQsuNZ9B1kB0bloVeO4pvVPOEiZb0yhUCo588Trd+03plVu6C7gG782U8ERdNMGr2oOHRG45pmDbnHW+OzOGYetlEVSe2IPnDLBlFb5t5+9CcNf34y+teWWNm0w7cxeKhYgmOltdxfm/rpcGq/Nme3o7P1vXL0k8AqITrYzPSlf2AOnTLCp625l0wab46/9w1cwXpue877mG4Ru8SybFRCbK1WuLxYgAoyO17BuaQFb71oz+/pWyiKpPTGA0xytzoCMklcAa2VXdq/jDw1XtItjATM3iCdePK7t/R84cdF18SxnpUpzfXv1ylXrFBDlE1MoNMt2d5u0amXTBi9uaZhSd8E1uHsd1636Jcj2a5Qf7IHTLLccchZ6gK1On3fjFojdasaXlYqu32qC7MmZ5m9JFC8GcJrVSg45LaLKE5vSM6ViAWMT+t43AGxcvcS1Msam+qXVyUHUvjxTKCLyVyLyvoj8uOmxnSJSEZGR+n/3RttMshG0BrohyMa/rZ7TL6/zRdUeU3rmi7cvRYdh38tSsYADJy66VsZ47TYfxuQgal82OfC/BvAFzeNPKaXuqP/3crjNIr/CyF/7zSHHnTP3Ol+U7dGV8T14Zxl7j1S0sySLhU7svG+N69om63ftB4A5xy0VC7ihe6Z0see6haFMDqL25ZlCUUodFJGV0TeFWhE0f92cTy11F9C82khjgozp/XHnzL3OF3V7dJOLdAOQnSKzgde0sQQwd0EqXaWKbuuwVittqL20kgP/qoj8FoDDAB5TSl3SvUhEtgHYBgDLly9HtVpt4ZTZNjY2Ftmxu6bGcetizRNT48Zr/tqpKgYPnkbX1DQ+tRgArgKdALpmnl+0YBL48LLx/UHOqWN7XdzO99C3X0HX1FXj83f/l79Bz3ULsXXdTdjw6XD2UzS2BwobVixCtVpF/+eWYfDgaXw0Oa17IYBJ/NnfDuNPXh6e1z7dddEdb9GCDmy7a1lufrei/D3KmqAB/M8BfBMz3bVvAvhjAF/RvVApNQhgEJjZ1Dhvm5E6RfX5P+zs1vbMyqWi8ZwDTx9DZVTBvFSUwsA/nEffr6wO7ZwmNq83nU8AvHu5Bq8lr969XMOxl8/gyQcWh9Ijt/n8fT09OPxebd4O8rbtc16Xvp4eoGvxnCqURzbnrwol73GkIVAduFLqPaXUlFJqGsBfAlgXbrPIryA10EFK2IKeM4zBRdN0ej8rv4Y54Gf7+Q+cuGjdRpv29a0t49COTfjnXb/KpWRzLlAPXESWKqUu1H/8dQA/dns9RS9IDbTfErag57QtfQuyTZrfBaiAazelVuupbT+/3wFGDkiSLc8ALiK7AXweQI+I/BTANwB8XkTuwEzn5wyA34mwjWTJbw2023ZcgN0sRptz2gwu2gZ521UKG8uvmgb8wqqntvn8fm80HJAkW54pFKXUVqXUUqVUQSm1XCn1HaXUl5VStymlPqOUuq+pN04Z4iyNu6G7MLv6nmm1uyCpEJvSt6D1zW5pDLfn4qyn9qr11rWPyAZnYuacn167W68VMKcSbErfgtY326QxnM812h7kfEE02vLYs0e1NeOdIphWal7bXztVxcDTxzhlnowYwMmaqdf6xIvH8WFt2piO0KVqnD1NryDvlq92uwk5n2vchEwUZtIyYQfLxrF018H0TWfw4Ol6lRCnzJMeAzhZM/VOdSvxNee4bXrJbkG+lXy1M/CPX5303PsyqmDpZ6B5YN9JdE1No7k0cqI2hceePYrte0bYIycADODkg9/BuOaA75WqcQtuuhmPtrNMnYHfVlQzSm1TVudHJ+qTq+ZqpGDYIyeAAZx8MPWSFy3omLMRQYPfagpTcAuaH291p/kky/mWlYrA1Ljra7K01C9Fgxs6kDXTvow771vT0kYKr52qula2BFklEWg9ACdZzte/ZRUWLfD+9WTNeL6xB06+uKUAgkyK8RqsGxqu4OcfTc57n80NwpTyKRY6cHVSaStC/By/VV4Ds/jwMgb+4TzOj06gQ0TbXtaM5xsDOIUi6EYKpsE63UbBDV6rJDboUj6FDsHktHvwbl5NMCo2A7MbPt0zuw6N8/UAa8aJAZwSZhqsOz86Ycxhdy9cYBVcdQOj41cnjftXAuayvjA097h1PWq3nHaU28VRdjGAU6JMg3XLSsVQNi9wfjO4ZcdLxteWIwyKzh606RuA22eLars4yi4OYlKidIN1jdSA7eCln+n9pmOWS8VIV/azrYhhTpv8YACnRPWtLWPbXZ+aV9nSmMHpVd3idxu1IMvuhsH2W8PG1UtaPlfc+5RScphCoUQ054PXLS2gf8v8QUnbdU78TPIxHROYmUIfVn7ZWWFS6i645t4bDpy46PvYzW3lrvX5wgBOsXMGmeqVq8Yg45X3DZInN62P4gx6h89+gAMnLgYqjXQer9AhKHQKalPuWzt49dS9AnQYe6NygDQ7mEKh2DS+2j+6ZyS0pVyDTvJpZgp6z7x+LtAO97rj1aYVPrZwwew65SZe7fZaBjfIDc1vGorSgwGcYtEcJEz8VJc0bgaV0Yl5O2H6zWmbzuvsK+tuMrp8s+l4YxM1HNqxyRjEBfBst1eADnJDi3NtdAoXAzjFwqYKw7bX7LwZNG/LbNqIIozzAnMDqKnnen2xoH1vqbtgvOkIgN/87E1Wu/u4PR7m3qicpp9+DOAUC69g4KfXrLsZKLiXArpVZpg2S9ZpDqCmnqsI5h2v0CkYm6jNuek0dIpAYWYA0ytt4RWgTevVeO2N6udxSg8OYlIs3Jai7bluIR6513tqfIPfHqPXwJ+uMmXj6iXYe6TiOnXddL7R8RqeeuiOOce79POPjAOYfpaItanMCWNvVE7TzwYGcIqFKUg8+cBt2LBiEXp6eqyPZbNFWzOvygxTBUbvzTe6Bkq3djiD6EqXGaCmdpmEPSOT0/SziwGcYuEWJKrVqq9j+e0xuvXYbXvnYbTDVhK5Z07TzyYGcIpNWEHCb4/RracctG7abztKxYJ20wtTe72wbpsABnDKKD83A7ee8vY9I9r32PaCbdux87416P/BUdSm3Sfy2PTgOduSGliFQm3PrTIjrgqMvrVlDHzpdpQ0JYZ+SyDjrNvmuirpxh44xW7e1//PLUOfj0HMIEw95TgrMBptaDX9EVfdNnv66ccATrHSBYXBg6eBrsWJBIUkKjBaHQvwW4UTVCvjAxQPBnCKlS4ofDQ5nWhQyFoFhulbw8bVS0JdUZEzNNOPAZxilYegEHWFiM3EozDSHXH19Ck4BnCKlDOYXW8opws7KCRVZhdX3tj5rWH9rv3adMdjzx7F9j0j2nXDva4PZ2imHwN4ymW53le7LnanoNAhc8rpFi3owCObwwsKSQ6+JZU3Nn2DaZ6m3//c0dnHba4PZ2imHwN4imW9CkC7LvaUwg3dBXQvXDAbFLbdtSzUz5Pk4JtpvZeoU0Q2O/7UphSeePE4uhcusL4+WRsfyBsG8BTLehWA22JPw1/fPPuz36n0Qc8bdRAdGq5AMH8dcSDavPHQcAVXPpy0eu2l8Zox0LfTOEReMICnWNYH/JIaBEvqvAP7TmqDt81GDa2e12uGp41lpWKmU3Z5xJmYKZb1dZqT2gE+bTvPK9invJpnPv7uM29ZzXwM44YuADauXsKt1TKGPfAUy3oVQFKDYEmd19Tz99oHs8G02bPX5spua63bamwokeWUXR55BnAR+SsAXwTwvlLql+qP3QhgD4CVAM4A+A2l1KXomplP7VAFkNQgWBLnbfWG67a5ciNBohvINp33wTvLs4HfK8FSLhUzn7LLI5se+F8D+FMA32t6bAeAv1NK7RKRHfWffz/85hGrAJJnmxdu9Ybrd3PlxnFtztvYi1OncZMZ2HdS+xrTHp+UPM8ArpQ6KCIrHQ/fD+Dz9T9/F8DfgwGcQhTnYJrbufyWcrZyw/WTCqmMTuCWHS/Naa/fzScA4IbuAn71M0uNwRsAfn51EkPDFXYkUihoDvwTSqkLAKCUuiAiHw+xTZRzcda/e50rylJO541Dtw+nm+aBxkZ7TUy99MNnP5iTotGpTSnmwVMq8kFMEdkGYBsALF++PPSa3ywZGxtLugmp5Lwuuw8ex7KuSaCr+dFJ7D54HBtWLAr13F7n6poax62LNW+cGm/p3/Jrp6oYPHgaXVPT+NTimeO9dfInePjOHgyfHUX1ylX8wr/oBGATzO2uzYYVi7DhP3xmThv+8fjZmfN7mRrHQ99+BVvX3YQNn4526V8v/D26JmgAf09EltZ730sBvG96oVJqEMAgAPT29io/m9e2o7x/fpPm6/LmhRrU7DYH15y+XAv9+nmd68PObmNlSSttGXj6GCqjCphzboWxE1dwaMfMJKdqtYr7nz5mlVYJcm0Gnj6GyuX5n93k3cs1HHv5DJ58IJmlf5vx92hG0DrwFwA8XP/zwwB+FE5ziOKtf/c6V1Q15bYVH7rz6yjA9445QapLotr5h4LxDOAishvAPwJYJSI/FZF/D2AXgHtE5BSAe+o/E4Uizok4Xudy246tFbY3Kef5S8UCCp36XrPfiTdBb4gsK0wPmyqUrYan7g65LUQA4q1/tzlXFKWcfmrGnedvDH7qUit+BlhNlSlesjITOA84E5NSKc769yRq7Vu5STXae8uOl7TVI7Y95OY2mPLszsW5sjQTOA8YwClT2mmxpTTsjdm80bJpNuffHL0wuwlHV4HLJ6UJ/zYoMxpBhostzQhzrMCU6++9+UZ8NDk9+7pL47VcX/O0YQ+cYtNq7znr66OHLeyxAl2u/T8+OwLnSrV5vuZpwwBOsXCb8Wg7OYeLLc0Xdv7ebYC0WZ6veZowhUKxcOs928r6+uhp15yi8lLqLsyuW+63/pzCwx44xSKM3nPW10dPO91N1uTKh5OzW7NVRiewfc8IHt0zgnJ9TRe39cspPAzgFIuwKiaAbK+P3qooq3D83EydW7g1r1f+/dfPzT6etY24s4YBnGIRVu85z+ujR71KYxg7++hw0DM6zIFTLKKakp4nYYwjuLFddyWIKG4MxB44xSjPvecwRF2F45aiMs36tNUp9qsekj32wIkyIo4qnL61ZRzasQlPPXQHAGD7nhGs37XfdVs1m9A8pVoJ/2TCAE6UEXGt0qib8frzq5ModMwN1cVCJ/7bQ3fgqYfumE2NmXraZZZ6RoIpFKKMiLoKx20ST21K4YbuAroXLtCe27SHKMBSzygxgBNlSFTjCLrA6zQ6XsM3fm3N7A2kMXjqXHoXyHepZ5wYwInIahJPqbtgVcbIwer4MAdORJ6VLMVCJ5RCpGWM5B974ETkOomnXE+DbN8zon2+1TLGqGaXttPa8SapD+B5+EsgSppppmzzZCvTAGcrZYxRzS6NetZqWqQ6hcIF/CnvhoYrsaz6ZzNTNooyRtPs0seePdrSZ4561mpapLoHzgX8Kc9eO1XF4y+fia0X6TX4GEWFiSn90pj4E/Qz52Xt+FQH8Lz8JRDp7H7zXGIdGFPqMuwKE5sFtIJ85jBWv8yCVKdQuIA/hcWUinBLUcSVvjCpXrmqfTzqDkycqUvbBbT8fua4Zq0mLdU9cC7gT2EwDWgdPvsB9h6paFMUABIfBOu5biHevVyb93jUHZg4U5fOtEyHiHbdFL+fOS8TilIdwPPyl0DRMgWk3W/8ZF6waB7oSnr8Zeu6m3CsKQcOxNOBiTt12ZyWCXMqfh4mFKU6gAP5+EugaHkNlNm+3uu5sG34dA+efGBx7B2YJPPH7LT5k/oATtQqU0Dq9Pi6noZBsCQ6MEmnLtlps5fqQUyiMJgGtLb+8grjQFdeBsF0srR7UtIDzUljD5zantvX8t6bb3T9ut7OX+XdZjlnoRecl9mWbkTFuFNGb2+vOnz4cGznS5tqtYqenp6km5E6vC56UV4X02BhWnvazRrXZf2u/do0V7lUxKEdmxJoWXRE5IhSqtf5OFMoRDnUDlPNOdGPKRSiXEpT8Au6YF1eZlu6YQ+cKIfSMsu5lVmfeR5obmAAJ8qhtAS/VlI5WaqWibjEKxUAAAWHSURBVApTKEQ5lJYJM62mcmyrZdp1XwEGcKI2ZBOwki4VHBquhLb2idd52rXcsKUUioicEZG3RWRERPJbH0iUIlnYCKXRRl3wDjuV0w4VNyZh5MA3KqXu0NUoElH8shCwdG0EZpY3CDuPnaaKm7BxEJOozWQhYJnaMq1U6GmNtFTcRKHVHLgC8IqIKAD/XSk16HyBiGwDsA0Ali9fjmq12uIps2tsbCzpJqQSr4te0OuybmlBuxlEz3ULU/P710ob/V6X/s8tw+DB0/hocnr2sUULOrDtrmWpuR5BtRrA1yulzovIxwG8KiInlFIHm19QD+qDwMxU+rxPmc775zfhddELcl223rVGO03+kXvXpOY6t9pGP5+jr6cH6Jq7LO8jm1mFAqXU+fr/3xeRHwJYB+Cg+7uIKEppKRF0E3cbk664iUrgAC4iHwPQoZS6XP/zZgB/GFrLiCiwLASsLLQx7VrpgX8CwA9FpHGc/6GU+l+htIqIiDwFDuBKqdMAbg+xLUREbSXqGaCciUlE5JNNYI5jBijrwImIfLCd6RrHhCoGcCIiH2wDcxwTqhjAiYh8sA3MccwAZQAnIvLBNjDHseY6AzgR5c7QcAXrd+3HLTtewvpd+32t1GgbmOPYcIJVKESUK61Wh/iZRRr1ZCUGcCLKFbdBSNtgm5ZZpEyhEFGuZGG5XVsM4ESUK+20PjgDOFEbaWVwLi/iqA6JC3PgRG2inTfvDVMWltu1xQBO1CbCGJzLi7QMQraKKRSiNtFOg3NkhwGcqE200+Ac2WEAJ2oT7TQ4R3aYAydqE+00OJdFUW/eoMMATtRG2mVwLmuSqgBiCoWIqEVxbN6gwx44EZFPznRJJaEKIAZwIiIfdOkSAaA0r426AogpFCIiH3TpEgVAHK+LowKIAZyIyAdTWkQBkW7eoMMUChGRD6acd7lUxKEdm2JtC3vgREQ+pGnCFHvgREQ+pGnCFAM4EZFPaZkwxRQKEVFGMYATEWUUAzgRUUYxgBMRZRQDOBFRRolSuhn8EZ1M5CKAs7GdMH16AFSTbkQK8bro8bro5fG63KyUWuJ8MNYAnnciclgp1Zt0O9KG10WP10WP1+UaplCIiDKKAZyIKKMYwOM1mHQDUorXRY/XRY/XpY45cCKijGIPnIgooxjAiYgyigE8ASLyn0REiUhP0m1JAxEZEJETInJMRH4oIqWk25QkEfmCiJwUkX8SkR1JtycNRGSFiBwQkXdE5LiI/F7SbUoDBvCYicgKAPcAOJd0W1LkVQC/pJT6DID/C+DxhNuTGBHpBPBnAP4VgF8EsFVEfjHZVqXCJIDHlFK/AOCzAH6X14UBPAlPAfjP0G9inUtKqVeUUpP1H18HsDzJ9iRsHYB/UkqdVkpdBfA/AdyfcJsSp5S6oJR6q/7nywDeAZD8gtwJYwCPkYjcB6CilDqadFtS7CsA/jbpRiSoDOAnTT//FAxUc4jISgBrAbyRbEuSxx15QiYi/xvAJzVP/QGArwHYHG+L0sHtuiilflR/zR9g5qvyM3G2LWVE8xi/rdWJyHUA9gJ4VCn1s6TbkzQG8JAppf6l7nERuQ3ALQCOiggwkyZ4S0TWKaX+X4xNTITpujSIyMMAvgjgbpXvyQk/BbCi6eflAM4n1JZUEZECZoL3M0qp55NuTxpwIk9CROQMgF6lVN5WVZtHRL4A4FsAPqeUuph0e5IkIgswM5B7N4AKgP8D4N8opY4n2rCEyUyv57sAPlBKPZp0e9KCOXBKgz8FsBjAqyIyIiJ/kXSDklIfzP0qgH2YGah7Nu/Bu249gC8D2FT/NzIiIvcm3aiksQdORJRR7IETEWUUAzgRUUYxgBMRZRQDOBFRRjGAExFlFAM4EVFGMYATEWXU/wc9Nl2dS6DdTwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(result.resid, result.fittedvalues)\n",
    "plt.grid(alpha = 0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we definitely see that the spread of the residuals increases with the fitted values.\n",
    "\n",
    "We have a case of heteroskedasticity\n",
    "\n",
    "Now heteroskedasticity does not mean our model should be thrown out\n",
    "just that there's more work that could be done to improve it.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# conclusion\n",
    "So as a conclusion, our newest model\n",
    "\n",
    "```Sales ~ sqrtTV + Radio + TV:Newspaper```\n",
    "\n",
    "is pretty good.\n",
    "\n",
    "We explain 93% of the variance in the Sales variable\n",
    "and besides a problem of heteroskedasticity\n",
    "all the assumptions of linear regression are basically satisfied.\n",
    "\n"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}