{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression diagnostics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example file shows how to use a few of the ``statsmodels`` regression diagnostic tests in a real-life context. You can learn about more tests and find out more information about the tests here on the [Regression Diagnostics page.](https://www.statsmodels.org/stable/diagnostic.html)\n",
    "\n",
    "Note that most of the tests described here only return a tuple of numbers, without any annotation. A full description of outputs is always included in the docstring and in the online ``statsmodels`` documentation. For presentation purposes, we use the ``zip(name,test)`` construct to pretty-print short descriptions in the examples below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimate a regression model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:34.726188Z",
     "iopub.status.busy": "2021-02-02T06:54:34.725507Z",
     "iopub.status.idle": "2021-02-02T06:54:34.959596Z",
     "shell.execute_reply": "2021-02-02T06:54:34.960084Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:34.965311Z",
     "iopub.status.busy": "2021-02-02T06:54:34.964696Z",
     "iopub.status.idle": "2021-02-02T06:54:35.702164Z",
     "shell.execute_reply": "2021-02-02T06:54:35.702573Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                Lottery   R-squared:                       0.348\n",
      "Model:                            OLS   Adj. R-squared:                  0.333\n",
      "Method:                 Least Squares   F-statistic:                     22.20\n",
      "Date:                Tue, 02 Feb 2021   Prob (F-statistic):           1.90e-08\n",
      "Time:                        06:54:35   Log-Likelihood:                -379.82\n",
      "No. Observations:                  86   AIC:                             765.6\n",
      "Df Residuals:                      83   BIC:                             773.0\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "===================================================================================\n",
      "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept         246.4341     35.233      6.995      0.000     176.358     316.510\n",
      "Literacy           -0.4889      0.128     -3.832      0.000      -0.743      -0.235\n",
      "np.log(Pop1831)   -31.3114      5.977     -5.239      0.000     -43.199     -19.424\n",
      "==============================================================================\n",
      "Omnibus:                        3.713   Durbin-Watson:                   2.019\n",
      "Prob(Omnibus):                  0.156   Jarque-Bera (JB):                3.394\n",
      "Skew:                          -0.487   Prob(JB):                        0.183\n",
      "Kurtosis:                       3.003   Cond. No.                         702.\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "from statsmodels.compat import lzip\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.stats.api as sms\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load data\n",
    "url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/HistData/Guerry.csv'\n",
    "dat = pd.read_csv(url)\n",
    "\n",
    "# Fit regression model (using the natural log of one of the regressors)\n",
    "results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()\n",
    "\n",
    "# Inspect the results\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normality of the residuals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jarque-Bera test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:35.709679Z",
     "iopub.status.busy": "2021-02-02T06:54:35.709271Z",
     "iopub.status.idle": "2021-02-02T06:54:35.713948Z",
     "shell.execute_reply": "2021-02-02T06:54:35.713573Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Jarque-Bera', 3.3936080248431666),\n",
       " ('Chi^2 two-tail prob.', 0.1832683123166337),\n",
       " ('Skew', -0.48658034311223375),\n",
       " ('Kurtosis', 3.003417757881633)]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['Jarque-Bera', 'Chi^2 two-tail prob.', 'Skew', 'Kurtosis']\n",
    "test = sms.jarque_bera(results.resid)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Omni test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:35.719914Z",
     "iopub.status.busy": "2021-02-02T06:54:35.719388Z",
     "iopub.status.idle": "2021-02-02T06:54:35.721894Z",
     "shell.execute_reply": "2021-02-02T06:54:35.722239Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Chi^2', 3.713437811597181), ('Two-tail probability', 0.15618424580304824)]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['Chi^2', 'Two-tail probability']\n",
    "test = sms.omni_normtest(results.resid)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Influence tests\n",
    "\n",
    "Once created, an object of class ``OLSInfluence`` holds attributes and methods that allow users to assess the influence of each observation. For example, we can compute and extract the first few rows of DFbetas by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:35.727115Z",
     "iopub.status.busy": "2021-02-02T06:54:35.725081Z",
     "iopub.status.idle": "2021-02-02T06:54:35.767952Z",
     "shell.execute_reply": "2021-02-02T06:54:35.767529Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.00301154,  0.00290872,  0.00118179],\n",
       "       [-0.06425662,  0.04043093,  0.06281609],\n",
       "       [ 0.01554894, -0.03556038, -0.00905336],\n",
       "       [ 0.17899858,  0.04098207, -0.18062352],\n",
       "       [ 0.29679073,  0.21249207, -0.3213655 ]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from statsmodels.stats.outliers_influence import OLSInfluence\n",
    "test_class = OLSInfluence(results)\n",
    "test_class.dfbetas[:5,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explore other options by typing ``dir(influence_test)``\n",
    "\n",
    "Useful information on leverage can also be plotted:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:35.772192Z",
     "iopub.status.busy": "2021-02-02T06:54:35.771427Z",
     "iopub.status.idle": "2021-02-02T06:54:35.938389Z",
     "shell.execute_reply": "2021-02-02T06:54:35.938756Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAGDCAYAAADHzQJ9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8oUlEQVR4nO3dfZzVdZ3//8fTEWRwkMEVEwdFMgRFTXI0W0xNbdFSIVpTyzTXvvxctbQLNikv8KLVsgsrlTK3VVaTFF3CdENNRfNqBccrIFYkTAY3cJMQGZKL1++Pz2fwMJyZOXPmnDkX87zfbufGOZ+r8zoXzOu8rxURmJmZWXXZrtQBmJmZWeE5wZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRVygjezTkkKSR9I7/9U0iUFvv4XJP2+kNds53k6jD3zdXbzeaZKuq2716kEPfXZWddtX+oArPJJWgZ8MSIeKnUslUjSXsAfgfsj4pMZ228DlkTE1BKFllVEnFPqGPJVybGbdZVL8FbxJNWUOoYCOUzS2O5eRFJV/3Cv9tdXSko4L1QJf5BWNJK2k3SRpFcl/Z+kOyXtnO77raTz2xz/gqSJ6f1Rkh6U9BdJiyV9JuO4WyRNk3S/pHeAj0n6pKQmSWskvS5paptrnyHptTSOSyQtk3RsZ3FmeU2LJJ2Q8Xh7SW9K+pCkfpJuS6+xWtKzkt7Xhbfsu8BVHbyf/0/SkvQ9mS1p94x9Iek8Sa8Ar0g6StJySf8iaaWkNyRNkPQJSf+TXuObGecfKumpNO43JF0vqW87cdwi6ar0/r2S1mbcNkv6Qrqvo8/w79LXsEbSfwN7d/C690pf39mS/gQ8nG7/p/TzeEvSHEnD0u2S9MP0df9V0ouS9m8be/p4cvp6V0j6pzbP+6ikL2Y83qoqWtKP0u/aGknzJX20nfhz/l5I+oakZklvp+/ZMen22jT2tyQtTONennHeVk0LbT6jQZJ+I2lVev5vJA1t8zq/LekJYB3w/kJ9dlZaTvBWTF8GJgBHArsDbwE3pPt+CZzWeqCk/YBhwH2SdgQeTI/ZNT3uRkmjM679WeDbwADg98A7wBlAPfBJ4J8lTci49o3A54AhwECgIcc427ojM25gHPBmRDwHnJleew/g74BzgJb2355t3ADso/SHRyZJRwNXA59JX8NrwIw2h00APgzslz7eDehH8lovBX4OnA4cDHwUuFTS+9NjNwFfAXYBPgIcA5zbWcARcWJE1EVEHfCPwP8Cv8vhM7wBWJ++ln9Kb505EtgXGJd+tt8EJgKDgcdJPhuAfwCOAPYh+T6cAvxf24tJOg74OvBxYASwzfveiWeBg4CdSV7nXZL6ZTkup++FpJHA+cAhETGA5Lu1LN19GUki3TvdfmYX4twO+HeS/197ps99fZtjPg9MIvn/tIrCf3ZWChHhm2/dupH8ETo2y/ZFwDEZj4cAG0j6fgwgScrD0n3fBn6R3j8FeLzNtX4GXJbevwWY3klM1wE/TO9fCtyRsa8/8G5rzB3FmeW6HwDeBvqnj28HLk3v/xPwJHBgF9+/vYBI35dzgafT7bcBU9P7/wZ8N+OcujTGvdLHARydsf8okj/kNenjAekxH844Zj4woZ2YLgT+M+NxAB/IeP+vanP8PsBK4KOdfYZATRr7qIx9/wr8vpP35/0Z2/4LODvj8XYkpc9hwNHA/wCHAdu1udaW2IFfANe0eQ2Zr/NRkr4lrfu/0F6M6f63gA+m96cCt3Xle5F+t1aS/NDo02bfUuC4jMeTgOXZPp/2PqOMfQcBb2U8fhS4IuNxwT4730p7cwneimkY8J9pteRqkkS6CXhfRLwN3Aecmh57KkmybD3vw63nped+jqRE2ur1zCeS9GFJj6TVkH8lKSXtku7ePfP4iFjH1iW6duNs+4IiYkm6/0RJ/YGTSEo6AP8BzAFmpFW+35XUJ4f3KdPPgfdJOrHN9t1JSu2tcaxNX0NmTcTrbc75v4jYlN5vLTH+OWN/C8kPBSTtk1bd/q+kNSR/tHchB5IGAr8GLomIx9PNHX2Gg0l+zGTG+xqdyzx+GPCjjGv/BRDQEBEPk5RQbwD+LOkmSTtlud5W34scY9hC0tfSJoK/pjEMJPt7ltP3Iv1uXUjy42ClpBl6rxkm71gl9Zf0MyVNVGuAx4B6bd13pe17W+jPzkrACd6K6XXg+Iioz7j1i4jmdP8dwGmSPgLUAo9knDe3zXl1EfHPGdduuwziL4HZwB4RMRD4KckffIA3gMw2x1qSqtJc42yrtZp+PLAw/cNMRGyIiMsjYj/g74ETSJoNchYRG4DLgSsz4gdYQfKHt/U17Ji+hswYu7M05DTgD8CIiNiJpPpbHZ+S9F8gee8fiYifZezq6DNcBWwkqbJutWcOMWa+vteB/6/N9Wsj4kmAiPhxRBwMjCYpmU/Ocr03OonhHZLanlZbfmCm7e3fIGkyGRQR9cBfyfKedeV7ERG/jIjDST7rAL6TY6zr2osV+BowkqT2ZieS5gvaxNr2vS30Z2cl4ARvhdIn7UzUetueJMl+W+91fhosaXzGOfeT/CG7AvhVRGxOt/+GpC3685L6pLdDJO3bwfMPAP4SEeslHUrSRt9qJkmJ+++VdBy7nK3/uHUWZ1szSNp5/5n3Su9I+pikA9KS0RqSqsxN2S/Rof8AdgCOy9j2S+AsSQdJ2oGkhP1MRCzL4/rZDCCJea2kUSSvLRffBnYELmizvd3PMK1VuAeYmpYu96NrbcqQfGZTWtuFJQ2UdHJ6/5C0RqcPSZJeT/bP4U7gC5L2S2tjLmuz/3lgYhrjB4CzM/YNIEl0q4DtJV0KZKslyPl7IWmkpKPTz3c9SQ1L63F3pq93kJIOcl/KEutnJdWkfQuObBNrC7BaSefRtq+zrWJ/dtZDnOCtUO4n+SPSepsK/IikVP2ApLeBp0k6gQEQEX8j+WNxLBmJMq2+/weSavsVJB23vkOS9NpzLnBF+jyXkvxBbL3eApI/iDNISkJvk7R1/i09pMM424qIN4CnSEpjv8rYtRvJj4k1JNX4c0na0VsnWPlpB/FnXn8TyR/hnTO2/Q64BLg7fQ17817zRiF8neRH0dskzQS/6vjwLU4jaet+S+/1pP9cDp/h+STNA/9L0l78710JNiL+M73ejLTa+WXg+HT3TulreIuk+vj/gO9lucZ/kfTVeBhYkv6b6YckfTX+DNzKe01IkFS5/xdJW/9rJAm5bRNJq3a/F23sAFwDvEnyvuxKUpMCyY/S10jmS3iA5EdgpguAE4HVJNXpszL2XUdSQ/YmyXf7t+3ECeT0/69bn531HEV0p1bPrPJIqiP5QzgiIv5Y4nDMukzSUSSd+IZ2cqj1Yi7BW68g6cS0SnFHktLcS7w3BMnMrOo4wVtvMZ6kunEFyZjnU8PVV2ZWxVxFb2ZmVoVcgjczM6tCTvBmZmZVqKpWZdpll11ir732KnUYZmZmPWL+/PlvRsTgbPuqKsHvtddezJs3r9RhmJmZ9QhJ7U4V7Cp6MzOzKuQEb2ZmVoWc4KvEK6+8Qr9+/Tj99NNLHYqZmZUBJ/gqcd5553HIIYeUOgwzMysTTvBVYMaMGdTX13PMMceUOhQzMysTTvAVbs2aNVx66aV8//vfL3UoZmZWRpzgK9wll1zC2WefzR577FHqUMzMrIxU1Tj43ub555/noYceoqmpqdShmJlZmXGCr2CPPvooy5YtY8899wRg7dq1bNq0iYULF/Lcc8+VODozMyulqlpNrrGxMXrTTHbr1q1jzZo1Wx5/73vfY9myZUybNo3Bg7POXGhmZlVE0vyIaMy2zyX4Cta/f3/69++/5XFdXR39+vVzcjczMyf4ajJ16tRSh2BmZmXCvejNzMyqkBO8mZlZFSpqgpd0nKTFkpZIuijL/lGSnpL0N0lfz7K/RlKTpN8UM04zM7NqU7QEL6kGuAE4HtgPOE3Sfm0O+wvwZeB77VzmAmBRsWI0MzOrVsUswR8KLImIpRHxLjADGJ95QESsjIhngQ1tT5Y0FPgkcHMRY6x4s5qaGXvNwwy/6D7GXvMws5qaSx2SmZmVgWIm+Abg9YzHy9NtuboO+Bdgc0cHSZokaZ6keatWrepykJVsVlMzU+55iebVLQTQvLqFKfe85CRvZmZFTfDKsi2nWXUknQCsjIj5nR0bETdFRGNENPa28d/XzllMy4ZNW21r2bCJa+csLlFEZmZWLoqZ4JcDmSugDAVW5HjuWOAkSctIqvaPlnRbYcOrfCtWt3Rpu5mZ9R7FTPDPAiMkDZfUFzgVmJ3LiRExJSKGRsRe6XkPR8TpxQu1Mu1eX9ul7WZm1nsULcFHxEbgfGAOSU/4OyNigaRzJJ0DIGk3ScuBrwIXS1ouaadixVRtJo8bSW2fmq221fapYfK4kSWKyMzMyoUXm6lws5qauXbOYlasbmH3+lomjxvJhDFd6ctoZmaVyovNVLEJYxqc0M3MbBueqtbMzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7wOaqrq9vqVlNTw5e+9KUt+3/3u98xatQo+vfvz8c+9jFee+21EkZrZma9nRN8jtauXbvl9uc//5na2lpOPvlkAN58800mTpzIlVdeyV/+8hcaGxs55ZRTShyxmZn1Zk7weZg5cya77rorH/3oRwG45557GD16NCeffDL9+vVj6tSpvPDCC/zhD38ocaRmZtZbOcHn4dZbb+WMM85ASpa8X7BgAR/84Ae37N9xxx3Ze++9WbBgQalCNDOzXs4Jvov+9Kc/MXfuXM4888wt29auXcvAgQO3Om7gwIG8/fbbPR2emZkZ4ATfZdOnT+fwww9n+PDhW7bV1dWxZs2arY5bs2YNAwYM6OnwzMzMACf4Lps+ffpWpXeA0aNH88ILL2x5/M477/Dqq68yevTong7PzMwMcILvkieffJLm5uYtvedbfepTn+Lll1/m7rvvZv369VxxxRUceOCBjBo1qkSRmplZb+cE3wW33norEydO3KbqffDgwdx9991861vfYtCgQTzzzDPMmDGjRFGamZmBIqLUMRRMY2NjzJs3r9RhmJmZ9QhJ8yOiMds+l+DNzMyqkBO8mZlZFXKCNzMzq0JO8GZmZlVo+1IHUMlmNTVz7ZzFrFjdwu71tUweN5IJYxpKHZaZmZkTfL5mNTUz5Z6XaNmwCYDm1S1MueclACd5MzMrOVfR5+naOYu3JPdWLRs2ce2cxSWKyMzM7D1O8HlasbqlS9vNzMx6khN8nnavr+3SdjMzs57kBJ+nyeNGUtunZqtttX1qmDxuZIkiMjMze4872eWptSOde9GbmVk5coLvhgljGpzQzcysLLmK3szMrAo5wZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqlBRE7yk4yQtlrRE0kVZ9o+S9JSkv0n6esb2PSQ9ImmRpAWSLihmnGZmZtWmaMPkJNUANwAfB5YDz0qaHRELMw77C/BlYEKb0zcCX4uI5yQNAOZLerDNuWZmZtaOYpbgDwWWRMTSiHgXmAGMzzwgIlZGxLPAhjbb34iI59L7bwOLAA84NzMzy1ExE3wD8HrG4+XkkaQl7QWMAZ4pTFhmZmbVr5gJXlm2RZcuINUBdwMXRsSado6ZJGmepHmrVq3KI0wzM7PqU8wEvxzYI+PxUGBFridL6kOS3G+PiHvaOy4iboqIxohoHDx4cN7BmpmZVZNiJvhngRGShkvqC5wKzM7lREkC/g1YFBE/KGKMZmZmValovegjYqOk84E5QA3wi4hYIOmcdP9PJe0GzAN2AjZLuhDYDzgQ+DzwkqTn00t+MyLuL1a8ZmZm1aSoq8mlCfn+Ntt+mnH/f0mq7tv6Pdnb8M3MzCwHnsnOzMysCjnBm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqkBO8mZlZFXKCNzMzq0JO8GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7wZmZmVcgJ3szMrAo5wefoqKOOol+/ftTV1VFXV8fIkSO37Fu3bh3nnnsuu+yyCwMHDuSII44oYaRmZmawfakDqCTXX389X/ziF7fZPmnSJDZu3MiiRYvYeeedef7553s+ODMzswxO8N20ePFiZs+ezfLly9lpp50AOPjgg0sclZmZ9Xauou+CKVOmsMsuuzB27FgeffRRAJ555hmGDRvGZZddxi677MIBBxzA3XffXdpAzcys13OCz9F3vvMdli5dSnNzM5MmTeLEE0/k1VdfZfny5bz88ssMHDiQFStWcP3113PmmWeyaNGiUodsZma9mBN8jj784Q8zYMAAdthhB84880zGjh3L/fffT21tLX369OHiiy+mb9++HHnkkXzsYx/jgQceKHXIZmbWiznB50kSEcGBBx5Y6lDMzMy24QSfg9WrVzNnzhzWr1/Pxo0buf3223nssccYN24cRxxxBHvuuSdXX301Gzdu5IknnuDRRx9l3LhxpQ7bzMx6Mfeiz8GGDRu4+OKL+cMf/kBNTQ2jRo1i1qxZW8bC//rXv+aLX/wi11xzDcOGDWP69OmMGjWqxFGbmVlvpogodQwF09jYGPPmzSt1GGZmZj1C0vyIaMy2z1X0ZmZmVcgJ3szMrAo5wZuZmVUhd7LL06ymZq6ds5gVq1vYvb6WyeNGMmFMQ6nDMjMzA5zg8zKrqZkp97xEy4ZNADSvbmHKPS8BOMmbmVlZcBV9Hq6ds3hLcm/VsmET185ZXKKIzMzMtuYEn4cVq1u6tN3MzKynOcHnYff62i5tNzMz62lO8HmYPG4ktX1qttpW26eGyeNGligiMzOzrRU1wUs6TtJiSUskXZRl/yhJT0n6m6Svd+XcUpowpoGrJx5AQ30tAhrqa7l64gHuYGdmZmWjaL3oJdUANwAfB5YDz0qaHRELMw77C/BlYEIe55bUhDENTuhmZla2ilmCPxRYEhFLI+JdYAYwPvOAiFgZEc8CG7p6rpmZmbWvmAm+AXg94/HydFtBz5U0SdI8SfNWrVqVV6BmZmbVppgJXlm25bp0Xc7nRsRNEdEYEY2DBw/OOTgzM7NqVswEvxzYI+PxUGBFD5xrZmbW6xUzwT8LjJA0XFJf4FRgdg+ca2Zm1usVrRd9RGyUdD4wB6gBfhERCySdk+7/qaTdgHnATsBmSRcC+0XEmmznFitWMzOzaqOIXJvFy19jY2PMmzev1GGYmZn1CEnzI6Ix2z7PZGdmZlaFnODNzMyqkBO8mZlZFXKCL4BXXnmFfv36cfrppwOwbNkyJFFXV7flduWVV5Y4SjMz602K1ou+NznvvPM45JBDttm+evVqtt/eb7GZmfU8l+C7acaMGdTX13PMMceUOhQzM7MtnOC7Yc2aNVx66aV8//vfz7p/2LBhDB06lLPOOos333yzh6MzM7PezAm+Gy655BLOPvts9thjj62277LLLjz77LO89tprzJ8/n7fffpvPfe5zJYrSzMx6IzcQ5+n555/noYceoqmpaZt9dXV1NDYm8w68733v4/rrr2fIkCGsWbOGnXbaqadDNTOzXsgJPk+PPvooy5YtY8899wRg7dq1bNq0iYULF/Lcc89tdayULI5XTbMGmplZefNUtXlat24da9as2fL4e9/7HsuWLWPatGksXbqU+vp6RowYwVtvvcW5557LypUreeSRR3okNjMz6x08VW0R9O/fn912223Lra6ujn79+jF48GCWLl3Kcccdx4ABA9h///3ZYYcduOOOO0odspmZ9SIuwZuZmVUol+DNzMx6GSd4MzOzKuQEb2ZmVoU8TK6AZjU1c+2cxaxY3cLu9bVMHjeSCWMaSh2WmZn1Qk7wBTKrqZkp97xEy4ZNADSvbmHKPS8BOMmbmVmPcxV9gVw7Z/GW5N6qZcMmrp2zuEQRmZlZb+YEXyArVrd0abuZmVkxOcEXyO71tV3abmZmVkxO8AUyedxIavvUbLWttk8Nk8eNLFFEZmbWm7mTXYG0dqRzL3ozMysHTvAFNGFMgxO6mZmVBVfRm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqUE4JXlJ/SZdI+nn6eISkE4obmpmZmeUr1xL8vwN/Az6SPl4OXFWUiMzMzKzbck3we0fEd4ENABHRAqhoUZmZmVm35Jrg35VUCwSApL1JSvRmZmZWhnKd6OYy4LfAHpJuB8YCXyhWUGZmZtY9OSX4iHhQ0nPAYSRV8xdExJtFjczMzMzyllOCl/Sh9O4b6b97ShoIvBYRG4sSmZmZmeUt1zb4G4GngZuAnwNPATOA/5H0D+2dJOk4SYslLZF0UZb9kvTjdP+LGT8kkPQVSQskvSzpDkn9uvTKzMzMerFcE/wyYExENEbEwcAY4GXgWOC72U6QVAPcABwP7AecJmm/NocdD4xIb5OAaem5DcCXgcaI2B+oAU7N/WWZmZn1brkm+FERsaD1QUQsJEn4Szs451BgSUQsjYh3SUr849scMx6YHomngXpJQ9J92wO1krYH+gMrcozVzMys18s1wS+WNE3SkentRpLq+R1Ix8Zn0QC8nvF4ebqt02Miohn4HvAnknb/v0bEAznGamZm1uvlmuC/ACwBLgS+AixNt20APtbOOdkmwolcjpE0iKR0PxzYHdhR0ulZn0SaJGmepHmrVq3q+FWYmZn1Ejkl+IhoiYjvR8SnImJCRHwvItZFxOaIWNvOacuBPTIeD2Xbavb2jjkW+GNErIqIDcA9wN+3E9tNad+AxsGDB+fycsrS6aefzpAhQ9hpp53YZ599uPnmm0sdkpmZVbBcF5sZIWmmpIWSlrbeOjntWWCEpOGS+pJ0kpvd5pjZwBlpb/rDSKri3yCpmj8sXeRGwDHAoi69sgozZcoUli1bxpo1a5g9ezYXX3wx8+fPL3VYZmZWobqy2Mw0YCNJlfx04D86OiEdH38+MIckOd8ZEQsknSPpnPSw+0mq+5eQDL87Nz33GWAm8BzwUhrnTbm/rMozevRodthhBwAkIYlXX321xFGZmVmlUkTbZvEsB0nzI+JgSS9FxAHptscj4qNFj7ALGhsbY968eaUOI2/nnnsut9xyCy0tLYwZM4bHHnuMurq6UodlZmZlKs3Pjdn25VqCXy9pO+AVSedL+hSwa8EiNABuvPFG3n77bR5//HEmTpy4pURvZmbWVbkm+AtJxqJ/GTgYOB04s0gx9Wo1NTUcfvjhLF++nGnTppU6HDMzy2LGjBnsu+++7Ljjjuy99948/vjjLFy4kMbGRgYNGsSgQYM49thjWbhwYcli7HQu+nRGus9ExGRgLXBW0aMyNm7c6DZ4M7My9OCDD/KNb3yDX/3qVxx66KG88UayTMuOO+7IzJkzGTZsGJs3b+aGG27g1FNP5cUXXyxJnJ2W4CNiE3Bw2pvdimDlypXMmDGDtWvXsmnTJubMmcMdd9zB0UcfXerQzMysjcsuu4xLL72Uww47jO22246GhgYaGhqor69nr732QhIRQU1NDUuWLClZnLmuB98E/FrSXcA7rRsj4p6iRNXLSGLatGmcc845bN68mWHDhnHdddcxfnzbmX3NzKyUNm3axLx58zjppJP4wAc+wPr165kwYQLXXnsttbW1ANTX17N27Vo2b97MFVdcUbJYc03wOwP/B2QWKYNkAhrrpsGDBzN37txSh2FmZp3485//zIYNG5g5cyaPP/44ffr0Yfz48Vx11VV8+9vfBmD16tW888473HrrrQwbNqxkseY0TK5SVPowOTMzK29vvfUWO++8M7fccgtnnpn0Nb/77ru56qqraGpq2urYzZs3M3jwYBYtWsSuuxZn4Fm3h8lJ2kfS7yS9nD4+UNLFhQzSzMys3A0aNIihQ4eSS7e0zZs3s27dOpqbm3sgsm3lOkzu58AU0pXjIuJFvD67mZn1QmeddRY/+clPWLlyJW+99RbXXXcdJ5xwAg8++CBNTU1s2rSJNWvW8NWvfpVBgwax7777liTOXBN8/4j47zbbNhY6mN5sVlMzY695mOEX3cfYax5mVlNpfvGZmVnHLrnkEg455BD22Wcf9t13X8aMGcO3vvUtVq9ezWmnncbAgQPZe++9WbJkCb/97W/p169fSeLMdara/yKZV/6uiPiQpH8Ezo6I44sdYFdUahv8rKZmptzzEi0bNm3ZVtunhqsnHsCEMQ0ljMzMzMpZIaaqPQ/4GTBKUjPJzHbndHiG5ezaOYu3Su4ALRs2ce2cxSWKyMzMKl2uw+Rei4hjJe0IbBcRbxczqN5mxeqWLm03MzPrTK4l+D9Kugk4jGS6Wiug3etru7TdzMysM7km+JHAQyRV9X+UdL2kw4sXVu8yedxIavvUbLWttk8Nk8eNLFFEZmZW6XKqoo+IFuBO4E5Jg4AfAXOBmg5PtJy0dqS7ds5iVqxuYff6WiaPG+kOdmZmFWBWU3NZ/v3OtQ0eSUcCpwDHA88CnylWUL3RhDENZfGFMDOz3LUdBdW8uoUp97wEUPK/6bnOZPdHkp7zjwP7R8RnIuLuYgZmZmZW7sp5FFSuJfgPRsSaokZiZmZWYcp5FFSunex281z0ZmZmWyvnUVCei97MzCxP5TwKKtcq+v4R8d9tVs/xXPRmZtarlfMoqFwT/JuS9gYCIJ2L/o2iRWUdKtchGWZmvVG5joLKNcGfB9zEe3PR/xH4XNGisnaV85AMMzMrHzm1wUfE0og4FhgMjIqIw4FPFTUyy6qch2SYmVn5yLWTHQAR8U7GQjNfLUI81on2hl40l8GQDDMzKx9dSvBtqPNDrNDaG3ohkup7MzMz6F6Cj4JFYTmbPG5k1l9WAa6mNzOzLTpM8JLelrQmy+1tYPceitEyTBjT0O4vq3KYOcnMzMpDh73oI2JATwViuWuor83a5l4OMyeZmVl56E4VvZVIOc+cZGZm5SHn5WKtfJTzzElmZlYenOArVLnOnGRmZuXBVfRmZmZVyCV461GeR9/MrGc4wVuP8Tz6ZmY9p6hV9JKOk7RY0hJJF2XZL0k/Tve/KOlDGfvqJc2U9AdJiyR9pJixWvF5Hn0zs55TtAQvqQa4ATge2A84TdJ+bQ47HhiR3iYB0zL2/Qj4bUSMAj4ILCpWrNYz2puIxxP0mJkVXjGr6A8FlkTEUgBJM4DxwMKMY8YD0yMigKfTUvsQ4B3gCOALABHxLvBuEWO1HrB7GU7Q4z4BZlatillF3wC8nvF4ebotl2PeD6wC/l1Sk6SbJe2Y7UkkTZI0T9K8VatWFS56K7hym6CntU9A8+oWgvf6BHjRHjOrBsVM8O2tiZLLMdsDHwKmRcQYkhL9Nm34ABFxU0Q0RkTj4MGDuxNvQcxqambsNQ8z/KL7GHvNw04WGSaMaeDqiQfQUF+LSKbcvXriASUrMbtPgJlVs2JW0S8H9sh4PBRYkeMxASyPiGfS7TNpJ8GXE/cS71w5TdDjPgFmVs2KWYJ/FhghabikvsCpwOw2x8wGzkh70x8G/DUi3oiI/wVel9Rad3sMW7fdlyWXCCtLe23/XrTHzKpB0RJ8RGwEzgfmkPSAvzMiFkg6R9I56WH3A0uBJcDPgXMzLvEl4HZJLwIHAf9arFgLxSXCylJufQLMzAqpqBPdRMT9JEk8c9tPM+4HcF475z4PNBYzvkIrx17i1j4v2mNm1cwz2RXQ5HEjt2qDB5cIy1059QkwMyskJ/gCconQzMzKhRN8gblEaGZm5cDLxZqZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpA72ZWAVzAzM7Nic4LvYZ6v3szMeoKr6HuY56s3M7Oe4ATfwzxfvZmZ9QQn+B7mFczMzKwnOMH3oFlNzax7d+M22z1fvZmZFZo72fWQtp3rWtXX9mHqSaPdwc7MzArKCb6AOhr+lq1zHcCOO2y/VXL3EDozMysEJ/gC6Wz4Wy6d6zyEzszMCsVt8AXS2fC3XDrXeQidmZkVihN8gXRWQp88biS1fWq22te2c52H0JmZWaE4wRdIZyX0CWMauHriATTU1yKgob6WqycesFXVu4fQmZlZobgNvkAmjxu5TS/5tiX0CWMaOmxLz+UaZmZmuXAJvkAmjGng0wc3UCMBUCPx6YM7TujZrtFZKd/MzCwXLsEXyKymZu6e38ymCAA2RXD3/GYah+3c5SRfSQndw/rMzMqTE3w3tSa45iwd4Vp7wJd7wss3SXtYn5lZ+XIVfTe0Jrhsyb1VufeAz3wNwXtJelZTc6fnelifmVn5cgm+G6bOXpB1drpMufaAL1VVd0dJurPn97A+M7Py5RJ8nmY1NbO6ZUOHx+TaA747peju6k6S9rA+M7Py5QSfp86qobvSA76UVd3dSdK5TN5jZmal4Sr6PHVUwr3ulIO6VL1eyqru7oy9z1xIx73ozczKixN8nnavr83auW5Q/z5dTnDtXasnqrq7m6QrbVifmVlv4QSfp/ZKvpedOLpg1+qpqm4naTOz6uMEn6dCVk/3ZFW3J6YxM+sdFOnMa9WgsbEx5s2bV+owylbbiWkgqSnwdLhmZpVJ0vyIaMy2zyX4AivnEnJ3xrybmVllcYIvoJ6aujXfHxGemMbMrPdwgu+Gtol23bsbi15C7s6PiFL21jczs55V1IluJB0nabGkJZIuyrJfkn6c7n9R0ofa7K+R1CTpN8WMMx/ZZp97a132me0KWULuzqQ4npjGzKz3KFqCl1QD3AAcD+wHnCZpvzaHHQ+MSG+TgGlt9l8ALCpWjN2RLdG2p5Al5O5Us3u9eTOz3qOYVfSHAksiYimApBnAeGBhxjHjgemRdOV/WlK9pCER8YakocAngW8DXy1inHnJtVReyBLyrKZmtpO2rDmfKdcfER7zbmbWOxQzwTcAr2c8Xg58OIdjGoA3gOuAfwEGFC/E/LXXnt1K6TGF6kXf2iSQLbm7mj035TzCwcys0IqZ4JVlW9vslPUYSScAKyNivqSjOnwSaRJJ9T577rlnHmHmZ/K4kXzlV89v84Igqfp+4qKjC/p87TUJ1EiuZs9BT41wMDMrF8XsZLcc2CPj8VBgRY7HjAVOkrQMmAEcLem2bE8SETdFRGNENA4ePLhQsXdqwpgGPnfYntv8QilWabq9JoHNEU5QOSjlin1mZqVQzAT/LDBC0nBJfYFTgdltjpkNnJH2pj8M+GtEvBERUyJiaETslZ73cEScXsRY83LVhAP44SkHddppbVZTM2OveZjhF93H2Gsezmudd6+93j2eA8DMepuiVdFHxEZJ5wNzgBrgFxGxQNI56f6fAvcDnwCWAOuAs4oVT6G1bc/9YTtLxBaqarjUC9JUOs8BYGa9jeeiz0NX5nQfe83DWRNLPu307iSWP8/Db2bVyHPRF1hX5nQvZNWwh7jlrydX7DMzKwdO8HnoStJ21XD58A8kM+tNijpVbbXqSoe33jQ9bCE6E5qZWWE4weehK0m7t0wPm21u/in3vOQkb2ZWIq6iz0NX23N7Q9VwJa41706LZlbNnODzlE/SruaEUmnjzD2znZlVO1fR95BZTc1MnvnCVlXYk2e+UDVV2JU2EY9ntjOzaucE30Muv3cBGzZtPefAhk3B5fcuKFFEhVVpnQkrrcbBzKyrXEXfQ95at6HT7ZVchV9p48w9fNHMqp0TfJmohjbhSupM6Kl/zazauYq+h9TX9ulwu9uEe1ZvGb5oZr2XS/A9ZOpJo5l81wts2PxeO3yf7cTUk0YDbhMuhUqqcTAz6yqX4PPU1VnbJoxp4NqTP7hVifHakz+4JcG01/Zb37+PZ4cz66a6urqtbjU1NXzpS1/asv/mm2/mAx/4AHV1dRx33HGsWLGihNGaFYZXk8tDMVYmy3bNPjWCYKtSv1dAM+ued955h/e9733cf//9HHHEEcydO5eTTz6ZRx55hBEjRnDBBRewcOFC5s6dW+pQzTrV0WpyLsHnoRjt5dnahHfsu/1Wyb0Qz2PW282cOZNdd92Vj370owDce++9nHzyyYwePZq+fftyySWX8Nhjj/Hqq6+WOFKz7nEbfB7aaxfPNuyqK9q2CQ+/6L6cnr+Sh9eZ9bRbb72VM844A0kARASZNZmt919++WX23nvvksRoVgguweehvfZyQUHbyHOZHa4Ui7x41TirVH/605+YO3cuZ5555pZtn/jEJ7jzzjt58cUXaWlp4YorrkAS69atK2GkZt3nBJ+H9sZKBxS0+jyX2eF6enidf1BYJZs+fTqHH344w4cP37LtmGOO4fLLL+fTn/40w4YNY6+99mLAgAEMHTq0hJGadZ8TfIF1NqytK8kql7HaPTm8blZTM1+784Wq/0Fh1Wv69Olbld5bnXfeebzyyiusXLmST3/602zcuJH999+/BBGaFY7b4PMwdXb788dvJzH8ovu2tIXDe9O31vfvw9r1G7d0nMtltrrOxmr31JSrrYl2UzujLoo1Xr8Sl6G18vTkk0/S3NzMySefvNX29evXs2TJEkaPHs3rr7/OpEmTuOCCCxg0aFCJIjUrDJfg87C6Jfu88gCbIt5bLe6uF7ZaQe6tdRsK3iu+pxZ5yZZoMxVrDvf2Oi56AiDrqltvvZWJEycyYMCArbavX7+ez372s9TV1XHooYfykY98hCuvvLJEUZoVjkvwRdQ2mbeneXULY695OK9e8Pku8tLVnvcdJdRizeE+q6kZkfRtaMuLwlhX/exnP8u6vb6+nhdffLGHozErPif4PAzq36fd1eHy1VpSzWeRma5OuZrPwjbtNQXUSEWbeOfaOYuzJnfRfkdHMzNLuIo+D5edOLqo1y/2ZDb59Lxvryng+5/5YNHawturNQgqZ4W9Qjr99NMZMmQIO+20E/vssw8333zzNsdcfvnlSOKhhx4qQYRmVk6c4PMwYUxDu6vDZeqznZLpZvPQ3UlzOpJPz/tSrL7WXjV8Qy+tnp8yZQrLli1jzZo1zJ49m4svvpj58+dv2f/qq68yc+ZMhgwZUsIozaxcOMHnaepJo+kodbcuJnPtP269wEyuapTfD4Nc5DKBTjYTxjTwxEVH88drPskTFx1d9FJ0T3UgrBSjR49mhx12AEASkraaTvX888/nO9/5Dn379i1ViBXF8ytYtXOCz9OEMQ1Z24chaSNuTYBtk2KuSb694WiFUCmJ02u2b+vcc8+lf//+jBo1iiFDhvCJT3wCgLvuuou+fftueWwd8/wK1hu4k103NOQxBn3yuJHbrBrX3rW7q72e8vn2vC8Fr9m+tRtvvJGf/OQnPPXUUzz66KPssMMOrF27lm9+85s88MADpQ6vYnh+BesNnOC7IVuy7qwk3Da5tp38Jpdr5KKznvJOnJWrpqaGww8/nNtuu41p06bx2muv8fnPf36r6VetYz05A6RZqTjBd0O+JeG2ybWrY9JzOd4llOq3ceNGXn31VebOncvy5cu58cYbAVi1ahWf+cxn+MY3vsE3vvGNEkdZnnpqBkizUnKC76ZClIS7co1cx7BXWgnFS952bOXKlTz88MOccMIJ1NbW8tBDD3HHHXfwy1/+kksvvZQNG96bl+GQQw7hBz/4Accff3wJIy5v+dS+mVUaJ/gKk2vJvBgllGIl4Xwm3ultJDFt2jTOOeccNm/ezLBhw7juuusYP378NsfW1NQwaNAg6urqShBpZaikfihm+XKCrzC5lswLXUIpZhJ2c0LnBg8ezNy5c3M6dtmyZcUNpkq4H4pVOw+TK5JijbHNdQx7oYeYFXPd+UprTjAzqwQuwRdBttLuhb96ngt/9Tz1tX2YetLovBNtV0rmhSyhFDMJu8OTmVnhuQRfBB0trbq6ZQOT73qhyyX61hqBr/zqeXbYfjsG9e/To5O/5Dv7XS4qZeIdM7NKUtQEL+k4SYslLZF0UZb9kvTjdP+Lkj6Ubt9D0iOSFklaIOmCYsZZaJ2Vajdsji5VbbeddWt1ywbWb9jMD085qEemjIXiJmHPWGdmVnhFq6KXVAPcAHwcWA48K2l2RCzMOOx4YER6+zAwLf13I/C1iHhO0gBgvqQH25xbttqrcs7U3o+AbD3Vy6ETWrF7HbvDU9d5aKGZdaSYbfCHAksiYimApBnAeCAzSY8HpkdEAE9Lqpc0JCLeAN4AiIi3JS0CGtqcW7ZymY42W9V2ez3V27tOT3dCcxIuHx5aaGadKWaCbwBez3i8nKR03tkxDaTJHUDSXsAY4JlsTyJpEjAJYM899+xuzF3WUSlq6uwFrG7ZsM05fbZT1qrt9krqNVLWxWfcCa33KodaHTMrb8VM8NnWO22bpTo8RlIdcDdwYUSsyfYkEXETcBNAY2Nj8ZZgyyKX+d5nNTVz+b0LeGtdkug76kXfXol8UwS1fWo865Zt4aGFZtaZYib45cAeGY+HAityPUZSH5LkfntE3FPEOPOWSymqK9Xa7bXd19f22XJtgEH9+3DZidl/JLhdtndo/a6smX8v77z8O95dtYwd9z2SD37umwC8++67fPazn2XevHm89tprPPLIIxx11FGlDdrMelQxe9E/C4yQNFxSX+BUYHabY2YDZ6S96Q8D/hoRb0gS8G/Aooj4QRFj7JZCl6Ky9VTvs514592NW1X1r9+wOev5XuO692j9rmxf93cM/Mgp1B3wcWraNP20rji32267lTBSMyuVopXgI2KjpPOBOUAN8IuIWCDpnHT/T4H7gU8AS4B1wFnp6WOBzwMvSXo+3fbNiLi/WPF21aym5qSBIUujQL5t49l6qq97d+OW6v1WLRs2MXX2gqKtIOdagPK35buyY19WrG5hh9XLGFH37pbtffv25cILLwSSuenNrPcp6kx2aUK+v822n2bcD+C8LOf9nuzt82VhVlMzX7vrBbL0e2u3A12u2lbpD7/ovqzHrW7ZwKym5oKvIFdNvbOr/YdK5nfl4oufYvny5SWOyMzKiWeyy8Pl9y5g0+bs/fnq+m1f0CTSUW1A28lyCjHbXDHnnO9Jbq4ws97OCT4PbavMM63uYF8+OqoNyLaCXHdnm6uW3tnV8kPFzCxfTvAFVuix6RPGNDCof5+cnqsQU74Wc875nlQtP1TMzPLlBJ+H1mFr2TSvbino8rAAl504OueS+YQxDTxx0dH88ZpP5jVPfbUs/FItP1Q6s3HjRtavX8+mTZvYtGkT69evZ+PGjQD87W9/Y/369UAybG79+vVEto4jZlaVnODzMPWk0fTZrv0+gNnae7uzPnxPLcbS2imtdfY8ivhcxdbVHyrd+XxK6aqrrqK2tpZrrrmG2267jdraWq666ioARo4cSW1tLc3NzYwbN47a2lpee+21EkdsZj1F1fSLvrGxMebNm9cjz3XxrJe4/ek/ZRslt0VDfS1PXHT0Nj3TIUk25ZQ4yzHG7vaCz/X8cnztZma5kDQ/Ihqz7SvqMLlq9sgfVnWY3OG99t5KmDe83GIsxHC9XGcRLLfXbmZWCE7wecqls1Zre285d/hqLeW2t7xtrjEWesx5Tybdcv58zMzy5Tb4PHXWWSuzvbdcO3xljhVvTy4xFmPMeU8m3XL9fMzMusMJPk/ZOnG1atsxLd+e6cXu+JWtlNzVGNu7TnfHnPdk0q2WkQNQuZ0FzazwnODzNGFMA58+uGGb+XRbE0NmNXI+veB7Yia2jkrDg/r3ybmTWTFK2z2ZdHtqlEKxefY+M8vkNvhuyNbRrr124q4sGws90wbd3vK00P6KdV25TndK29kW3inmXPJd/XzKkTsLmlkmJ/huKGY7cU+0QU8eN3Kb4WGtupIYsl2nEKXtaki6PcmdBc0sk6vou6GY7cQ90QbdWjXdnlwTQ7VUcVc6dxY0s0xO8N1QzHbinmqDnjCmgYYCJIbuTpFr3VdNnQXNrPuc4LuhmCXXniwVOzFUB9ekmFkmT1VrQOEnqjEzs+LzVLXWKXdoMzOrLq6iNzMzq0JO8GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7wZmZmVcgJ3szMrAo5wZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRVSRJQ6hoKRtAp4rYefdhfgzR5+zmrj97D7/B52n9/DwvD72H1deQ+HRcTgbDuqKsGXgqR5EdFY6jgqmd/D7vN72H1+DwvD72P3Feo9dBW9mZlZFXKCNzMzq0JO8N13U6kDqAJ+D7vP72H3+T0sDL+P3VeQ99Bt8GZmZlXIJXgzM7Mq5ASfJ0nHSVosaYmki0odT6WRtIekRyQtkrRA0gWljqlSSaqR1CTpN6WOpVJJqpc0U9If0u/kR0odU6WR9JX0//LLku6Q1K/UMVUCSb+QtFLSyxnbdpb0oKRX0n8H5XNtJ/g8SKoBbgCOB/YDTpO0X2mjqjgbga9FxL7AYcB5fg/zdgGwqNRBVLgfAb+NiFHAB/H72SWSGoAvA40RsT9QA5xa2qgqxi3AcW22XQT8LiJGAL9LH3eZE3x+DgWWRMTSiHgXmAGML3FMFSUi3oiI59L7b5P8QW0obVSVR9JQ4JPAzaWOpVJJ2gk4Avg3gIh4NyJWlzSoyrQ9UCtpe6A/sKLE8VSEiHgM+EubzeOBW9P7twIT8rm2E3x+GoDXMx4vx8kpb5L2AsYAz5Q4lEp0HfAvwOYSx1HJ3g+sAv49beq4WdKOpQ6qkkREM/A94E/AG8BfI+KB0kZV0d4XEW9AUhgCds3nIk7w+VGWbR6OkAdJdcDdwIURsabU8VQSSScAKyNifqljqXDbAx8CpkXEGOAd8qwS7a3SNuLxwHBgd2BHSaeXNipzgs/PcmCPjMdDcXVUl0nqQ5Lcb4+Ie0odTwUaC5wkaRlJM9HRkm4rbUgVaTmwPCJaa5BmkiR8y92xwB8jYlVEbADuAf6+xDFVsj9LGgKQ/rsyn4s4wefnWWCEpOGS+pJ0Jpld4pgqiiSRtHkuiogflDqeShQRUyJiaETsRfIdfDgiXGrqooj4X+B1SSPTTccAC0sYUiX6E3CYpP7p/+1jcEfF7pgNnJnePxP4dT4X2b5g4fQiEbFR0vnAHJLeor+IiAUlDqvSjAU+D7wk6fl02zcj4v7ShWS92JeA29Mf7EuBs0ocT0WJiGckzQSeIxkh04RntMuJpDuAo4BdJC0HLgOuAe6UdDbJj6eT87q2Z7IzMzOrPq6iNzMzq0JO8GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb1YEkkLS9zMef13S1B6O4VFJjen9+yXVd/N6RxVjxTpJT7az/RZJ/5jnNTuMNR2rTetnolSWbQdJeipdJe1FSafkE49ZKXgcvFlx/A2YKOnqiHizqydL2j4iNhYqmIj4RKGu1ZF84o6IUsx49hVJa0imVP02MBfYP8u2ZcAZEfGKpN2B+ZLmeDEaqwQuwZsVx0aSiT6+0naHpGGSfpeWCH8nac90+y2SfiDpEeA76eNpkh6RtFTSkena0Ysk3ZJxvWmS5qWlzMuzBSNpmaRdJJ0j6fn09sf0uZD0D2lJ9TlJd6VrBCDpOCVrpP8emNjOtb+QnnMv8ICkHdM4n00XbxmfHjda0n+nz/2ipBHp9rXpv5J0vaSFku4jY4GN1vjT+42SHk3vHyrpyfR5nsyYjS4zviMzXnOTpAHp7Im7kCxx+tuIeKCdbf8TEa8ARMQKkilDB2d7H8zKjRO8WfHcAHxO0sA2268HpkfEgcDtwI8z9u0DHBsRX0sfDwKOJvmhcC/wQ2A0cICkg9JjvhURjcCBwJGSDmwvoIj4aUQcBBxCMgf7D9LEeXH6vB8C5gFfldQP+DlwIvBRYLcOXutHgDMj4mjgWyTT5h4CfAy4VsnqbOcAP0qfvzF9/kyfAkYCBwD/j9zmMv8DcES6SMylwL9mOebrwHnp834UaJF0IfAmyXt/nKSPZ9uWeRFJhwJ9gVdziMus5FxFb1YkEbFG0nSSEmFLxq6P8F5p+D+A72bsuysiNmU8vjciQtJLwJ8j4iUASQuAvYDngc9ImkTy/3kIsB/wYifh/YgkCd+rZFW6/YAn0mbovsBTwCiSBUReSZ/zNmBSO9d7MCJa17T+B5JFcL6ePu4H7Jle81tK1rC/p/W6GY4A7khf/wpJD3fyGgAGAremtQEB9MlyzBMkP2RuT593uaQfpe/r1IiYmra/P5RlG+lrH0LyWZ0ZEV6a1yqCS/BmxXUdcDbQ0frimfNFv9Nm39/Sfzdn3G99vL2k4SQl1GPSGoH7SBJquyR9ARgGtFbniyRBH5Te9ouIs7PE1pHMuAV8OuN6e0bEooj4JXASyY+dOZKOznKd9p5vI+/9vcp8fVcCj0TE/iQ1Ddu89oi4BvgiUAs8LWlUpHN0R8TU9N/Itg1A0k4k7+vFEfF0x2+DWflwgjcrorRUeydJkm/1JMnqbwCfA37fjafYiSS5/lXS+4DjOzpY0sEkPwhOzyiJPg2MlfSB9Jj+kvYhqf4eLmnv9LjTcoxpDvCl1hKwpDHpv+8HlkbEj0lWy2rblPAYcKqkmrTE/LGMfcuAg9P7n87YPhBoTu9/oZ3XvHdEvBQR3yFpfhiV4+tAyeIz/0nSpHJXrueZlQMneLPi+z5J561WXwbOkvQiyYp6F+R74Yh4gWTlrgXAL0iqoztyPrAz8Eja6ezmiFhFkhzvSGN6GhgVEetJquTvSzvZvZZjWFeSVJW/KOnl9DHAKcDLSlYPHAVMb3PefwKvAC8B00h6sbe6HPiRpMeBzCaM7wJXS3qCZGXHbC6U9LKkF0hqD/4rx9cB8BmSpoMvZHTUO6gL55uVjFeTMzMzq0IuwZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRVygjczM6tCTvBmZmZVyAnezMysCv3/HOzPxZ3N6fgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from statsmodels.graphics.regressionplots import plot_leverage_resid2\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "fig = plot_leverage_resid2(results, ax = ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other plotting options can be found on the [Graphics page.](https://www.statsmodels.org/stable/graphics.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multicollinearity\n",
    "\n",
    "Condition number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:35.944346Z",
     "iopub.status.busy": "2021-02-02T06:54:35.943769Z",
     "iopub.status.idle": "2021-02-02T06:54:35.946519Z",
     "shell.execute_reply": "2021-02-02T06:54:35.946861Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "702.1792145490062"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(results.model.exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heteroskedasticity tests\n",
    "\n",
    "Breush-Pagan test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:35.951851Z",
     "iopub.status.busy": "2021-02-02T06:54:35.951072Z",
     "iopub.status.idle": "2021-02-02T06:54:35.955846Z",
     "shell.execute_reply": "2021-02-02T06:54:35.956170Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Lagrange multiplier statistic', 4.893213374093957),\n",
       " ('p-value', 0.08658690502352209),\n",
       " ('f-value', 2.503715946256434),\n",
       " ('f p-value', 0.08794028782673029)]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['Lagrange multiplier statistic', 'p-value',\n",
    "        'f-value', 'f p-value']\n",
    "test = sms.het_breuschpagan(results.resid, results.model.exog)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goldfeld-Quandt test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:35.960963Z",
     "iopub.status.busy": "2021-02-02T06:54:35.960184Z",
     "iopub.status.idle": "2021-02-02T06:54:35.964915Z",
     "shell.execute_reply": "2021-02-02T06:54:35.965254Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('F statistic', 1.1002422436378152), ('p-value', 0.3820295068692507)]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['F statistic', 'p-value']\n",
    "test = sms.het_goldfeldquandt(results.resid, results.model.exog)\n",
    "lzip(name, test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linearity\n",
    "\n",
    "Harvey-Collier multiplier test for Null hypothesis that the linear specification is correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:54:35.969519Z",
     "iopub.status.busy": "2021-02-02T06:54:35.968732Z",
     "iopub.status.idle": "2021-02-02T06:54:35.974569Z",
     "shell.execute_reply": "2021-02-02T06:54:35.974213Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('t value', -1.0796490077823977), ('p value', 0.28346392475408044)]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name = ['t value', 'p value']\n",
    "test = sms.linear_harvey_collier(results)\n",
    "lzip(name, test)"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}