{
 "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "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:                Fri, 10 Jul 2020   Prob (F-statistic):           1.90e-08\n",
      "Time:                        05:46:36   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",
      "Warnings:\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",
    "import statsmodels.datasets\n",
    "dat = statsmodels.datasets.get_rdataset(\"Guerry\", \"HistData\", cache=True).data\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": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Jarque-Bera', 3.3936080248431755),\n",
       " ('Chi^2 two-tail prob.', 0.18326831231663288),\n",
       " ('Skew', -0.48658034311223436),\n",
       " ('Kurtosis', 3.0034177578816346)]"
      ]
     },
     "execution_count": 3,
     "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": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Chi^2', 3.7134378115971938), ('Two-tail probability', 0.15618424580304724)]"
      ]
     },
     "execution_count": 4,
     "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": {
    "collapsed": true
   },
   "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": 5,
     "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": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAGDCAYAAADHzQJ9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfZzVdZ3//8fTEWRwkMEVEwdFMgRFTXI0W0xNbdFSIVpTyzTXvvxctbQLNikv8KLVsgsrlTK3VVaTFF3CdENNRfNqBccrIFYkSAY3cJMQGZKL1++Pz2fwMJ6ZOXPmnDkX87zfbufGOZ+r8zoXzOu8rxURmJmZWXXZrtQBmJmZWeE5wZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRVygjezTkkKSR9I7/9U0iUFvv4XJP2+kNds53k6jD3zdXbzeaZKuq2716kEPfXZWddtX+oArPJJWgZ8MSIeKnUslUjSXsAfgfsj4pMZ228DlkTE1NJEll1EnFPqGPJVybGbdZVL8FbxJNWUOoYCOUzS2O5eRFJV/3Cv9tdXSko4L1QJf5BWNJK2k3SRpFcl/Z+kOyXtnO77raTz2xz/gqSJ6f1Rkh6U9BdJiyV9JuO4WyRNk3S/pLeBj0n6pKQmSWslvSZpaptrnyFpeRrHJZKWSTq2szizvKZFkk7IeLy9pDckfUhSP0m3pddYI+lZSe/rwlv2XeCqDt7P/ydpSfqezJa0e8a+kHSepFeAVyQdJWmFpH+RtErS65ImSPqEpP9Jr/HNjPMPlfRUGvfrkq6X1LedOG6RdFV6/15J6zJuWyR9Id3X0Wf4d+lrWCvpv4G9O3jde6Wv72xJfwIeTrf/U/p5vClpjqRh6XZJ+mH6uv8q6UVJ+7eNPX08OX29KyX9U5vnfVTSFzMeb1MVLelH6XdtraT5kj7aTvw5fy8kfUNSs6S30vfsmHR7bRr7m5IWpnGvyDhvm6aFNp/RIEm/kbQ6Pf83koa2eZ3flvQEsB54f6E+OystJ3grpi8DE4Ajgd2BN4Eb0n2/BE5rPVDSfsAw4D5JOwIPpsfsmh53o6TRGdf+LPBtYADwe+Bt4AygHvgk8M+SJmRc+0bgc8AQYCDQkGOcbd2RGTcwDngjIp4DzkyvvQfwd8A5QEvHb9E2bgD2UfrDI5Oko4Grgc+kr2E5MKPNYROADwP7pY93A/qRvNZLgZ8DpwMHAx8FLpX0/vTYzcBXgF2AjwDHAOd2FnBEnBgRdRFRB/wj8L/A73L4DG8ANqSv5Z/SW2eOBPYFxqWf7TeBicBg4HGSzwbgH4AjgH1Ivg+nAP/X9mKSjgO+DnwcGAG8533vxLPAQcDO6eu8S1K/LMfl9L2QNBI4HzgkIgaQfLeWpbsvI0mke6fbz+xCnNsB/07y/2vP9Lmvb3PM54FJJP+fVlP4z85KISJ8861bN5I/Qsdm2b4IOCbj8RBgI0nfjwEkSXlYuu/bwC/S+6cAj7e51s+Ay9L7twDTO4npOuCH6f1LgTsy9vUH3mmNuaM4s1z3A8BbQP/08e3Apen9fwKeBA7s4vu3FxDp+3Iu8HS6/TZganr/34DvZpxTl8a4V/o4gKMz9h9F8oe8Jn08ID3mwxnHzAcmtBPThcB/ZjwO4AMZ7/9VbY7fB1gFfLSzzxCoSWMflbHvX4Hfd/L+vD9j238BZ2c83o6k9DkMOBr4H+AwYLs219oaO/AL4Jo2ryHzdT5K0rekdf8X2osx3f8m8MH0/lTgtq58L9Lv1iqSHxp92uxbChyX8XgSsCLb59PeZ5Sx7yDgzYzHjwJXZDwu2GfnW2lvLsFbMQ0D/jOtllxDkkg3A++LiLeA+4BT02NPJUmWred9uPW89NzPkZRIW72W+USSPizpkbQa8q8kpaRd0t27Zx4fEevZtkTXbpxtX1BELEn3nyipP3ASSUkH4D+AOcCMtMr3u5L65PZWbfVz4H2STmyzfXeSUntrHOvS15BZE/Fam3P+LyI2p/dbS4x/ztjfQvJDAUn7pFW3/ytpLckf7V3IgaSBwK+BSyLi8XRzR5/hYJIfM5nxLqdzmccPA36Uce2/AAIaIuJhkhLqDcCfJd0kaacs19vme5FjDFtJ+lraRPDXNIaBZH/PcvpepN+tC0l+HKySNEPvNsPkHauk/pJ+pqSJai3wGFCvbfuutH1vC/3ZWQk4wVsxvQYcHxH1Gbd+EdGc7r8DOE3SR4Ba4JGM8+a2Oa8uIv4549ptl0H8JTAb2CMiBgI/JfmDD/A6kNnmWEtSVZprnG21VtOPBxamf5iJiI0RcXlE7Af8PXACSbNBziJiI3A5cGVG/AArSf7wtr6GHdPXkBljd5aGnAb8ARgRETuRVH+r41OS/gsk7/0jEfGzjF0dfYargU0kVdat9swhxszX9xrw/7W5fm1EPAkQET+OiIOB0SQl88lZrvd6JzG8TVLb02rrD8y0vf0bJE0mgyKiHvgrWd6zrnwvIuKXEXE4yWcdwHdyjHV9e7ECXwNGktTe7ETSfEGbWNu+t4X+7KwEnOCtUPqknYlab9uTJNlv693OT4Mljc84536SP2RXAL+KiC3p9t+QtEV/XlKf9HaIpH07eP4BwF8iYoOkQ0na6FvNJClx/72SjmOXs+0ft87ibGsGSTvvP/Nu6R1JH5N0QFoyWktSlbk5+yU69B/ADsBxGdt+CZwl6SBJO5CUsJ+JiGV5XD+bASQxr5M0iuS15eLbwI7ABW22t/sZprUK9wBT09LlfnStTRmSz2xKa7uwpIGSTk7vH5LW6PQhSdIbyP453Al8QdJ+aW3MZW32Pw9MTGP8AHB2xr4BJIluNbC9pEuBbLUEOX8vJI2UdHT6+W4gqWFpPe7O9PUOUtJB7ktZYv2spJq0b8GRbWJtAdYo6Tza9nW2VezPznqIE7wVyv0kf0Rab1OBH5GUqh+Q9BbwNEknMAAi4m8kfyyOJSNRptX3/0BSbb+SpOPWd0iSXnvOBa5In+dSkj+IrddbQPIHcQZJSegtkrbOv6WHdBhnWxHxOvAUSWnsVxm7diP5MbGWpBp/Lkk7eusEKz/tIP7M628m+SO8c8a23wGXAHenr2Fv3m3eKISvk/woeoukmeBXHR++1Wkkbd1v6t2e9J/L4TM8n6R54H9J2ov/vSvBRsR/ptebkVY7vwwcn+7eKX0Nb5JUH/8f8L0s1/gvkr4aDwNL0n8z/ZCkr8afgVt5twkJkir3/yJp619OkpDbNpG0avd70cYOwDXAGyTvy64kNSmQ/ChdTjJfwgMkPwIzXQCcCLRWp8/K2HcdSQ3ZGyTf7d+2EyeQ0/+/bn121nMU0Z1aPbPKI6mO5A/hiIj4Y6njMesqSUeRdOIb2tmx1nu5BG+9gqQT0yrFHUlKcy/x7hAkM7Oq4wRvvcV4kurGlSRjnk8NV1+ZWRVzFb2ZmVkVcgnezMysCjnBm5mZVaGqWpVpl112ib322qvUYZiZmfWI+fPnvxERg7Ptq6oEv9deezFv3rxSh2FmZtYjJLU7VbCr6M3MzKqQE7yZmVkVcoKvEq+88gr9+vXj9NNPL3UoZmZWBpzgq8R5553HIYccUuowzMysTDjBV4EZM2ZQX1/PMcccU+pQzMysTDjBV7i1a9dy6aWX8v3vf7/UoZiZWRlxgq9wl1xyCWeffTZ77LFHqUMxM7MyUlXj4Hub559/noceeoimpqZSh2JmZmXGCb6CPfrooyxbtow999wTgHXr1rF582YWLlzIc889V+LozMyslKpqNbnGxsboTTPZrV+/nrVr1259/L3vfY9ly5Yxbdo0Bg/OOnOhmZlVEUnzI6Ix2z6X4CtY//796d+//9bHdXV19OvXz8ndzMyc4KvJ1KlTSx2CmZmVCfeiNzMzq0JO8GZmZlWoqAle0nGSFktaIumiLPtHSXpK0t8kfT3L/hpJTZJ+U8w4zczMqk3RErykGuAG4HhgP+A0Sfu1OewvwJeB77VzmQuARcWK0czMrFoVswR/KLAkIpZGxDvADGB85gERsSoingU2tj1Z0lDgk8DNRYyx4s1qambsNQ8z/KL7GHvNw8xqai51SGZmVgaKmeAbgNcyHq9It+XqOuBfgC0dHSRpkqR5kuatXr2661FWsFlNzUy55yWa17QQQPOaFqbc85KTvJmZFTXBK8u2nGbVkXQCsCoi5nd2bETcFBGNEdHY28Z/XztnMS0bN2+zrWXjZq6ds7hEEZmZWbkoZoJfAWSugDIUWJnjuWOBkyQtI6naP1rSbYUNr/KtXNPSpe1mZtZ7FDPBPwuMkDRcUl/gVGB2LidGxJSIGBoRe6XnPRwRpxcv1Mq0e31tl7abmVnvUbQEHxGbgPOBOSQ94e+MiAWSzpF0DoCk3SStAL4KXCxphaSdihVTtZk8biS1fWq22Vbbp4bJ40aWKCIzMysXXmymws1qaubaOYtZuaaF3etrmTxuJBPGdKUvo5mZVSovNlPFJoxpcEI3M7P38FS1ZmZmVcgJ3szMrAo5wZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRVygs9RXV3dNreamhq+9KUvbd3/u9/9jlGjRtG/f38+9rGPsXz58hJGa2ZmvZ0TfI7WrVu39fbnP/+Z2tpaTj75ZADeeOMNJk6cyJVXXslf/vIXGhsbOeWUU0ocsZmZ9WZO8HmYOXMmu+66Kx/96EcBuOeeexg9ejQnn3wy/fr1Y+rUqbzwwgv84Q9/KHGkZmbWWznB5+HWW2/ljDPOQEqWvF+wYAEf/OAHt+7fcccd2XvvvVmwYEGpQjQzs17OCb6L/vSnPzF37lzOPPPMrdvWrVvHwIEDtzlu4MCBvPXWWz0dnpmZGeAE32XTp0/n8MMPZ/jw4Vu31dXVsXbt2m2OW7t2LQMGDOjp8MzMzAAn+C6bPn36NqV3gNGjR/PCCy9sffz222/z6quvMnr06J4Oz8zMDHCC75Inn3yS5ubmrb3nW33qU5/i5Zdf5u6772bDhg1cccUVHHjggYwaNapEkZqZWW/nBN8Ft956KxMnTnxP1fvgwYO5++67+da3vsWgQYN45plnmDFjRomiNDMzA0VEqWMomMbGxpg3b16pwzAzM+sRkuZHRGO2fS7Bm5mZVSEneDMzsyrkBG9mZlaFnODNzMyq0PalDqCSzWpq5to5i1m5poXd62uZPG4kE8Y0lDosMzMzJ/h8zWpqZso9L9GycTMAzWtamHLPSwBO8mZmVnKuos/TtXMWb03urVo2bubaOYtLFJGZmdm7nODztHJNS5e2m5mZ9SQn+DztXl/bpe1mZmY9yQk+T5PHjaS2T80222r71DB53MgSRWRmZvYud7LLU2tHOveiNzOzcuQE3w0TxjQ4oZuZWVlyFb2ZmVkVcoI3MzOrQk7wZmZmVcgJ3szMrAo5wZuZmVWhoiZ4ScdJWixpiaSLsuwfJekpSX+T9PWM7XtIekTSIkkLJF1QzDjNzMyqTdGGyUmqAW4APg6sAJ6VNDsiFmYc9hfgy8CENqdvAr4WEc9JGgDMl/Rgm3PNzMysHcUswR8KLImIpRHxDjADGJ95QESsiohngY1ttr8eEc+l998CFgEecG5mZpajYib4BuC1jMcryCNJS9oLGAM8U5CozMzMeoFiJnhl2RZduoBUB9wNXBgRa9s5ZpKkeZLmrV69Oo8wzczMqk8xE/wKYI+Mx0OBlbmeLKkPSXK/PSLuae+4iLgpIhojonHw4MF5B2tmZlZNipngnwVGSBouqS9wKjA7lxMlCfg3YFFE/KCIMZqZmVWlovWij4hNks4H5gA1wC8iYoGkc9L9P5W0GzAP2AnYIulCYD/gQODzwEuSnk8v+c2IuL9Y8ZqZmVWToq4mlybk+9ts+2nG/f8lqbpv6/dkb8M3MzOzHHgmOzMzsyrkBG9mZlaFnODNzMyqkBO8mZlZFXKCNzMzq0JO8GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7wZmZmVcgJ3szMrAo5wZuZmVUhJ3gzM7Mq5ASfo6OOOop+/fpRV1dHXV0dI0eO3Lpv/fr1nHvuueyyyy4MHDiQI444ooSRmpmZwfalDqCSXH/99Xzxi198z/ZJkyaxadMmFi1axM4778zzzz9fgujMzMze5QTfTYsXL2b27NmsWLGCnXbaCYCDDz64xFGZmVlv5yr6LpgyZQq77LILY8eO5dFHHwXgmWeeYdiwYVx22WXssssuHHDAAdx9992lDdTMzHo9J/gcfec732Hp0qU0NzczadIkTjzxRF599VVWrFjByy+/zMCBA1m5ciXXX389Z555JosWLSp1yGZm1os5wefowx/+MAMGDGCHHXbgzDPPZOzYsdx///3U1tbSp08fLr74Yvr27cuRRx7Jxz72MR544IFSh2xmZr2YE3yeJBERHHjggaUOxczM7D2c4HOwZs0a5syZw4YNG9i0aRO33347jz32GOPGjeOII45gzz335Oqrr2bTpk088cQTPProo4wbN67UYZuZWS/mXvQ52LhxIxdffDF/+MMfqKmpYdSoUcyaNWvrWPhf//rXfPGLX+Saa65h2LBhTJ8+nVGjRpU4ajMz680UEaWOoWAaGxtj3rx5pQ7DzMysR0iaHxGN2fa5it7MzKwKOcGbmZlVISd4MzOzKuROdnma1dTMtXMWs3JNC7vX1zJ53EgmjGkodVhmZmaAE3xeZjU1M+Wel2jZuBmA5jUtTLnnJQAneTMzKwuuos/DtXMWb03urVo2bubaOYtLFJGZmdm2nODzsHJNS5e2m5mZ9TQn+DzsXl/bpe1mZmY9zQk+D5PHjaS2T80222r71DB53MgSRWRmZratoiZ4ScdJWixpiaSLsuwfJekpSX+T9PWunFtKE8Y0cPXEA2ior0VAQ30tV088wB3szMysbBStF72kGuAG4OPACuBZSbMjYmHGYX8BvgxMyOPckpowpsEJ3czMylYxS/CHAksiYmlEvAPMAMZnHhARqyLiWWBjV881MzOz9hUzwTcAr2U8XpFuK+i5kiZJmidp3urVq/MK1MzMrNoUM8Ery7Zcl67L+dyIuCkiGiOicfDgwTkHZ2ZmVs2KmeBXAHtkPB4KrOyBc83MzHq9Yib4Z4ERkoZL6gucCszugXPNzMx6vaL1oo+ITZLOB+YANcAvImKBpHPS/T+VtBswD9gJ2CLpQmC/iFib7dxixWpmZlZtFJFrs3j5a2xsjHnz5pU6DDMzsx4haX5ENGbb55nszMzMqpATvJmZWRVygjczM6tCTvAF8Morr9CvXz9OP/10AJYtW4Yk6urqtt6uvPLKEkdpZma9SdF60fcm5513Hocccsh7tq9Zs4btt/dbbGZmPc8l+G6aMWMG9fX1HHPMMaUOxczMbCsn+G5Yu3Ytl156Kd///vez7h82bBhDhw7lrLPO4o033ujh6MzMrDdzgu+GSy65hLPPPps99thjm+277LILzz77LMuXL2f+/Pm89dZbfO5znytRlGZm1hu5gThPzz//PA899BBNTU3v2VdXV0djYzLvwPve9z6uv/56hgwZwtq1a9lpp516OlQzM+uFnODz9Oijj7Js2TL23HNPANatW8fmzZtZuHAhzz333DbHSsnieNU0a6CZmZU3T1Wbp/Xr17N27dqtj7/3ve+xbNkypk2bxtKlS6mvr2fEiBG8+eabnHvuuaxatYpHHnmkR2IzM7PewVPVFkH//v3Zbbfdtt7q6uro168fgwcPZunSpRx33HEMGDCA/fffnx122IE77rij1CGbmVkv4hK8mZlZhXIJ3szMrJdxgjczM6tCTvBmZmZVyMPkCmhWUzPXzlnMyjUt7F5fy+RxI5kwpqHUYZmZWS/kBF8gs5qamXLPS7Rs3AxA85oWptzzEoCTvJmZ9ThX0RfItXMWb03urVo2bubaOYtLFJGZmfVmTvAFsnJNS5e2m5mZFZMTfIHsXl/bpe1mZmbF5ARfIJPHjaS2T80222r71DB53MgSRWRmZr2ZO9kVSGtHOveiNzOzcuAEX0ATxjQ4oZuZWVlwFb2ZmVkVcoI3MzOrQk7wZmZmVcgJ3szMrArllOAl9Zd0iaSfp49HSDqhuKGZmZlZvnItwf878DfgI+njFcBVRYnIzMzMui3XBL93RHwX2AgQES2AihaVmZmZdUuuCf4dSbVAAEjam6REb2ZmZmUo14luLgN+C+wh6XZgLPCFYgVlZmZm3ZNTgo+IByU9BxxGUjV/QUS8UdTIzMzMLG85JXhJH0rvvp7+u6ekgcDyiNhUlMjMzMwsb7m2wd8IPA3cBPwceAqYAfyPpH9o7yRJx0laLGmJpIuy7JekH6f7X8z4IYGkr0haIOllSXdI6telV2ZmZtaL5ZrglwFjIqIxIg4GxgAvA8cC3812gqQa4AbgeGA/4DRJ+7U57HhgRHqbBExLz20Avgw0RsT+QA1wau4vy8zMrHfLNcGPiogFrQ8iYiFJwl/awTmHAksiYmlEvENS4h/f5pjxwPRIPA3USxqS7tseqJW0PdAfWJljrGZmZr1ergl+saRpko5MbzeSVM/vQDo2PosG4LWMxyvSbZ0eExHNwPeAP5G0+/81Ih7IMVYzM7NeL9cE/wVgCXAh8BVgabptI/Cxds7JNhFO5HKMpEEkpfvhwO7AjpJOz/ok0iRJ8yTNW716dScvw8zMrHfIKcFHREtEfD8iPhUREyLiexGxPiK2RMS6dk5bAeyR8Xgo761mb++YY4E/RsTqiNgI3AP8fTux3ZT2DWgcPHhwLi+nLJ1++ukMGTKEnXbaiX322Yebb7651CGZmVkFy3WxmRGSZkpaKGlp662T054FRkgaLqkvSSe52W2OmQ2ckfamP4ykKv51kqr5w9JFbgQcAyzq0iurMFOmTGHZsmWsXbuW2bNnc/HFFzN//vxSh2VmZhWqK4vNTAM2kVTJTwf+o6MT0vHx5wNzSJLznRGxQNI5ks5JD7ufpLp/Ccnwu3PTc58BZgLPAS+lcd6U+8uqPKNHj2aHHXYAQBKSePXVV0sclZmZVSpFtG0Wz3KQND8iDpb0UkQckG57PCI+WvQIu6CxsTHmzZtX6jDydu6553LLLbfQ0tLCmDFjeOyxx6irqyt1WGZmVqbS/NyYbV+uJfgNkrYDXpF0vqRPAbsWLEID4MYbb+Stt97i8ccfZ+LEiVtL9GZmZl2Va4K/kGQs+peBg4HTgTOLFVRvVlNTw+GHH86KFSuYNm1aqcMxM7MsZsyYwb777suOO+7I3nvvzeOPP87ChQtpbGxk0KBBDBo0iGOPPZaFCxeWLMZO56JPZ6T7TERMBtYBZxU9KmPTpk1ugzczK0MPPvgg3/jGN/jVr37FoYceyuuvJ8u07LjjjsycOZNhw4axZcsWbrjhBk499VRefPHFksTZaQk+IjYDB6e92a0IVq1axYwZM1i3bh2bN29mzpw53HHHHRx99NGlDs3MzNq47LLLuPTSSznssMPYbrvtaGhooKGhgfr6evbaay8kERHU1NSwZMmSksWZ63rwTcCvJd0FvN26MSLuKUpUvYwkpk2bxjnnnMOWLVsYNmwY1113HePHt53Z18zMSmnz5s3MmzePk046iQ984ANs2LCBCRMmcO2111JbWwtAfX0969atY8uWLVxxxRUlizXXBL8z8H9AZpEySCagsW4aPHgwc+fOLXUYZmbWiT//+c9s3LiRmTNn8vjjj9OnTx/Gjx/PVVddxbe//W0A1qxZw9tvv82tt97KsGHDShZrTsPkKkWlD5MzM7Py9uabb7Lzzjtzyy23cOaZSV/zu+++m6uuuoqmpqZtjt2yZQuDBw9m0aJF7LprcQaedXuYnKR9JP1O0svp4wMlXVzIIM3MzMrdoEGDGDp0KLl0S9uyZQvr16+nubm5ByJ7r1yHyf0cmEK6clxEvIjXZzczs17orLPO4ic/+QmrVq3izTff5LrrruOEE07gwQcfpKmpic2bN7N27Vq++tWvMmjQIPbdd9+SxJlrgu8fEf/dZtumQgfTm81qambsNQ8z/KL7GHvNw8xqKs0vPjMz69gll1zCIYccwj777MO+++7LmDFj+Na3vsWaNWs47bTTGDhwIHvvvTdLlizht7/9Lf369StJnLlOVftfJPPK3xURH5L0j8DZEXF8sQPsikptg5/V1MyUe16iZePmrdtq+9Rw9cQDmDCmoYSRmZlZOSvEVLXnAT8DRklqJpnZ7pyOT7FcXTtn8TbJHaBl42aunbO4RBGZmVmly3WY3PKIOFbSjsB2EfFWMYPqbVauaenSdjMzs87kWoL/o6SbgMNIpqu1Atq9vrZL283MzDqTa4IfCTxEUlX/R0nXSzq8eGH1LpPHjaS2T80222r71DB53MgSRWRmZpUupyr6iGgB7gTulDQI+BEwF6jp8ETLSWtHumvnLGblmhZ2r69l8riR7mBnZlYBZjU1l+Xf71zb4JF0JHAKcDzwLPCZYgXVG00Y01AWXwgzM8td21FQzWtamHLPSwAl/5ue60x2fyTpOf84sH9EfCYi7i5qZGZmZmWunEdB5VqC/2BErC1qJGZmZhWmnEdB5drJbjfPRW9mZratch4F5bnozczM8lTOo6ByraLvHxH/3Wb1HM9Fb2ZmvVo5j4LKNcG/IWlvIADSuehfL1pU1qFyHZJhZtYblesoqFwT/HnATbw7F/0fgc8VLSprVzkPyTAzs/KRUxt8RCyNiGOBwcCoiDgc+FRRI7OsynlIhpmZlY9cO9kBEBFvZyw089UixGOdaG/oRXMZDMkwM7Py0aUE34Y6P8QKrb2hFyKpvjczM4PuJfgoWBSWs8njRmb9ZRXganozM9uqwwQv6S1Ja7Pc3gJ276EYLcOEMQ3t/rIqh5mTzMysPHTYiz4iBvRUIJa7hvrarG3u5TBzkpmZlYfuVNFbiZTzzElmZlYecl4u1spHOc+cZGZm5cEJvkKV68xJZmZWHlxFb2ZmVoVcgrce5Xn0zcx6hhO89RjPo29m1nOKWkUv6ThJiyUtkXRRlv2S9ON0/4uSPpSxr17STEl/kLRI0keKGasVn+fRNzPrOUVL8MMi1OcAABygSURBVJJqgBuA44H9gNMk7dfmsOOBEeltEjAtY9+PgN9GxCjgg8CiYsVqPaO9iXg8QY+ZWeEVs4r+UGBJRCwFkDQDGA8szDhmPDA9IgJ4Oi21DwHeBo4AvgAQEe8A7xQxVusBu5fhBD3uE2Bm1aqYVfQNwGsZj1ek23I55v3AauDfJTVJulnSjtmeRNIkSfMkzVu9enXhoreCK7cJelr7BDSvaSF4t0+AF+0xs2pQzATf3poouRyzPfAhYFpEjCEp0b+nDR8gIm6KiMaIaBw8eHB34i2IWU3NjL3mYYZfdB9jr3nYySLDhDENXD3xABrqaxHJlLtXTzygZCVm9wkws2pWzCr6FcAeGY+HAitzPCaAFRHxTLp9Ju0k+HLiXuKdK6cJetwnwMyqWTFL8M8CIyQNl9QXOBWY3eaY2cAZaW/6w4C/RsTrEfG/wGuSWutuj2Hbtvuy5BJhZWmv7d+L9phZNShago+ITcD5wBySHvB3RsQCSedIOic97H5gKbAE+DlwbsYlvgTcLulF4CDgX4sVa6G4RFhZyq1PgJlZIRV1opuIuJ8kiWdu+2nG/QDOa+fc54HGYsZXaOXYS9za50V7zKyaeSa7Apo8buQ2bfDgEmG5K6c+AWZmheQEX0AuEZqZWblwgi8wlwjNzKwceLlYMzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcie7EvAKZmZmVmxO8D3M89WbmVlPcBV9D/N89WZm1hOc4HuY56s3M7Oe4ATfw7yCmZmZ9QQn+B40q6mZ9e9ses92z1dvZmaF5k52PaRt57pW9bV9mHrSaHewMzOzgnKCL6COhr9l61wHsOMO22+T3D2EzszMCsEJvkA6G/6WS+c6D6EzM7NCcRt8gXQ2/C2XznUeQmdmZoXiBF8gnZXQJ48bSW2fmm32te1c5yF0ZmZWKE7wBdJZCX3CmAaunngADfW1CGior+XqiQdsU/XuIXRmZlYoboMvkMnjRr6nl3zbEvqEMQ0dtqXncg0zM7NcuARfIBPGNPDpgxuokQCokfj0wR0n9GzX6KyUb2ZmlguX4AtkVlMzd89vZnMEAJsjuHt+M43Ddu5ykq+khO5hfWZm5ckJvptaE1xzlo5wrT3gyz3h5ZukPazPzKx8uYq+G1oTXLbk3qrce8Bnvobg3SQ9q6m503M9rM/MrHy5BN8NU2cvyDo7XaZce8CXqqq7oyTd2fN7WJ+ZWflyCT5Ps5qaWdOyscNjcu0B351SdHd1J0l7WJ+ZWflygs9TZ9XQXekBX8qq7u4k6Vwm7zEzs9JwFX2eOirhXnfKQV2qXi9lVXd3xt5nLqTjXvRmZuXFCT5Pu9fXZu1cN6h/ny4nuPau1RNV3d1N0pU2rM/MrLdwgs9TeyXfy04cXbBr9VRVt5O0mVn1cYLPUyGrp3uyqtsT05iZ9Q6KdOa1atDY2Bjz5s0rdRhlq+3ENJDUFHg6XDOzyiRpfkQ0ZtvnEnyBlXMJuTtj3s3MrLI4wRdQT03dmu+PCE9MY2bWezjBd0PbRLv+nU1FLyF350dEKXvrm5lZzyrqRDeSjpO0WNISSRdl2S9JP073vyjpQ23210hqkvSbYsaZj2yzz725PvvMdoUsIXdnUhxPTGNm1nsULcFLqgFuAI4H9gNOk7Rfm8OOB0akt0nAtDb7LwAWFSvG7siWaNtTyBJyd6rZvd68mVnvUcwq+kOBJRGxFEDSDGA8sDDjmPHA9Ei68j8tqV7SkIh4XdJQ4JPAt4GvFjHOvORaKi9kCXlWUzPbSVvXnM+U648Ij3k3M+sdipngG4DXMh6vAD6cwzENwOvAdcC/AAOKGGPe2mvPbqX0mEL1om9tEsiW3F3NnptyHuFgZlZoxUzwyrKtbXbKeoykE4BVETFf0lEdPok0iaR6nz333DOfOPMyedxIvvKr59/zgiCp+n7ioqML+nztNQnUSK5mz0FPjXAwMysXxexktwLYI+PxUGBljseMBU6StAyYARwt6bZsTxIRN0VEY0Q0Dh48uFCxd2rCmAY+d9ie7/mFUqzSdHtNAlsinKByUMoV+8zMSqGYCf5ZYISk4ZL6AqcCs9scMxs4I+1Nfxjw14h4PSKmRMTQiNgrPe/hiDi9iLHm5aoJB/DDUw7qtNParKZmxl7zMMMvuo+x1zyc1zrvXnu9ezwHgJn1NkWroo+ITZLOB+YANcAvImKBpHPS/T8F7gc+ASwB1gNnFSueQmvbnvvDdpaILVTVcKkXpKl0ngPAzHobz0Wfh67M6T72moezJpZ82undSSx/noffzKqR56IvsK7M6V7IqmEPcctfT67YZ2ZWDpzg89CVpO2q4fLhH0hm1psUdaraatWVDm+9aXrYQnQmNDOzwnCCz0NXknZvmR4229z8U+55yUnezKxEXEWfh6625/aGquFKXGvenRbNrJo5wecpn6RdzQml0saZe2Y7M6t2rqLvIbOampk884VtqrAnz3yhaqqwK20iHs9sZ2bVzgm+h1x+7wI2bt52zoGNm4PL711QoogKq9I6E1ZajYOZWVe5ir6HvLl+Y6fbK7kKv9LGmXv4oplVOyf4MlENbcKV1JnQU/+aWbVzFX0Pqa/t0+F2twn3rN4yfNHMei+X4HvI1JNGM/muF9i45d12+D7biaknjQbcJlwKlVTjYGbWVS7B56mrs7ZNGNPAtSd/cJsS47Unf3Brgmmv7be+fx/PDmfWTXV1ddvcampq+NKXvrR1/80338wHPvAB6urqOO6441i5cmUJozUrDK8ml4dirEyW7Zp9agTBNqV+r4Bm1j1vv/0273vf+7j//vs54ogjmDt3LieffDKPPPIII0aM4IILLmDhwoXMnTu31KGadaqj1eRcgs9DMdrLs7UJ79h3+22SeyGex6y3mzlzJrvuuisf/ehHAbj33ns5+eSTGT16NH379uWSSy7hscce49VXXy1xpGbd4zb4PLTXLp5t2FVXtG0THn7RfTk9fyUPrzPrabfeeitnnHEGkgCICDJrMlvvv/zyy+y9994lidGsEFyCz0N77eWCgraR5zI7XCkWefGqcVap/vSnPzF37lzOPPPMrds+8YlPcOedd/Liiy/S0tLCFVdcgSTWr19fwkjNus8JPg/tjZUOKGj1eS6zw/X08Dr/oLBKNn36dA4//HCGDx++ddsxxxzD5Zdfzqc//WmGDRvGXnvtxYABAxg6dGgJIzXrPif4AutsWFtXklUuY7V7cnjdrKZmvnbnC1X/g8Kq1/Tp07cpvbc677zzeOWVV1i1ahWf/vSn2bRpE/vvv38JIjQrHLfB52Hq7Pbnj99OYvhF921tC4d3p2+t79+HdRs2be04l8tsdZ2N1e6pKVdbE+3mdkZdFGu8fiUuQ2vl6cknn6S5uZmTTz55m+0bNmxgyZIljB49mtdee41JkyZxwQUXMGjQoBJFalYYLsHnYU1L9nnlATZHvLta3F0vbLOC3JvrNxa8V3xPLfKSLdFmKtYc7u11XPQEQNZVt956KxMnTmTAgAHbbN+wYQOf/exnqaur49BDD+UjH/kIV155ZYmiNCscl+CLqG0yb0/zmhbGXvNwXr3g813kpas97ztKqMWaw31WUzMi6dvQlheFsa762c9+lnV7fX09L774Yg9HY1Z8TvB5GNS/T7urw+WrtaSazyIzXZ1yNZ+FbdprCqiRijbxzrVzFmdN7qL9jo5mZpZwFX0eLjtxdFGvX+zJbPLped9eU8D3P/PBorWFt1drEFTOCnuFdPrppzNkyBB22mkn9tlnH26++eb3HHP55ZcjiYceeqgEEZpZOXGCz8OEMQ3trg6Xqc92SqabzUN3J83pSD4970ux+lp71fANvbR6fsqUKSxbtoy1a9cye/ZsLr74YubPn791/6uvvsrMmTMZMmRICaM0s3LhBJ+nqSeNpqPU3bqYzLX/uO0CM7mqUX4/DHKRywQ62UwY08ATFx3NH6/5JE9cdHTRS9E91YGwUowePZoddtgBAElI2mY61fPPP5/vfOc79O3bt1QhVhTPr2DVzgk+TxPGNGRtH4akjbg1AbZNirkm+faGoxVCpSROr9n+Xueeey79+/dn1KhRDBkyhE984hMA3HXXXfTt23frY+uY51ew3sCd7LqhIY8x6JPHjXzPqnHtXbu72uspn2/P+1Lwmu3buvHGG/nJT37CU089xaOPPsoOO+zAunXr+OY3v8kDDzxQ6vAqhudXsN7ACb4bsiXrzkrCbZNr28lvcrlGLjrrKe/EWblqamo4/PDDue2225g2bRrLly/n85///DbTr1rHenIGSLNScYLvhnxLwm2Ta1fHpOdyvEso1W/Tpk28+uqrzJ07lxUrVnDjjTcCsHr1aj7zmc/wjW98g2984xsljrI89dQMkGal5ATfTYUoCXflGrmOYa+0EoqXvO3YqlWrePjhhznhhBOora3loYce4o477uCXv/wll156KRs3vjsvwyGHHMIPfvADjj/++BJGXN7yqX0zqzRO8BUm15J5MUooxUrC+Uy809tIYtq0aZxzzjls2bKFYcOGcd111zF+/Pj3HFtTU8OgQYOoq6srQaSVoZL6oZjlywm+wuRaMi90CaWYSdjNCZ0bPHgwc+fOzenYZcuWFTeYKuF+KFbtPEyuSIo1xjbXMeyFHmJWzHXnK605wcysErgEXwTZSrsX/up5LvzV89TX9mHqSaPzTrRdKZkXsoRSzCTsDk9mZoXnEnwRdLS06pqWjUy+64Uul+hbawS+8qvn2WH77RjUv0+PTv6S7+x3uaiUiXfMzCpJURO8pOMkLZa0RNJFWfZL0o/T/S9K+lC6fQ9Jj0haJGmBpAuKGWehdVaq3bglulS13XbWrTUtG9mwcQs/POWgHpkyFoqbhD1jnZlZ4RWtil5SDXAD8HFgBfCspNkRsTDjsOOBEentw8C09N9NwNci4jlJA4D5kh5sc27Zaq/KOVN7PwKy9VQvh05oxe517A5PXeehhWbWkWK2wR8KLImIpQCSZgDjgcwkPR6YHhEBPC2pXtKQiHgdeB0gIt6StAhoaHNu2cplOtpsVdvt9VRv7zo93QnNSbh8eGihmXWmmAm+AXgt4/EKktJ5Z8c0kCZ3AEl7AWOAZ7I9iaRJwCSAPffcs5shd11HpaipsxewpmXje87ps52yVm23V1KvkbIuPuNOaL1XOdTqmFl5K2aCz7beadss1eExkuqAu4ELI2JttieJiJuAmwAaGxuLtwRbFrnM9z6rqZnL713Am+uTRN9RL/r2SuSbI6jtU+NZt2wrDy00s84UM8GvAPbIeDwUWJnrMZL6kCT32yPiniLGmbdcSlFdqdZur+2+vrbP1msDDOrfh8tOzP4jwe2yvUPrd2Xt/Ht5++Xf8c7qZey475F88HPfBOCdd97hs5/9LPPmzWP58uU88sgjHHXUUaUN2sx6VDF70T8LjJA0XFJf4FRgdptjZgNnpL3pDwP+GhGvSxLwb8CiiPhBEWPslkKXorL1VO+znXj7nU3bVPVv2Lgl6/le47r3aP2ubF/3dwz8yCnUHfBxato0/bSuOLfbbruVMFIzK5WileAjYpOk84E5QA3wi4hYIOmcdP9PgfuBTwBLgPXAWenpY4HPAy9Jej7d9s2IuL9Y8XbVrKbmpIEhS6NAvm3j2Xqqr39n09bq/VYtGzczdfaCoq0g51qA8rf1u7JjX1auaWGHNcsYUffO1u19+/blwgsvBJK56c2s9ynqTHZpQr6/zbafZtwP4Lws5/2e7O3zZWFWUzNfu+sFsvR7a7cDXa7aVukPv+i+rMetadnIrKbmgq8gV029s6v9h0rmd+Xii59ixYoVJY7IzMqJZ7LLw+X3LmDzluz9+er6bV/QJNJRbUDbyXIKMdtcMeec70lurjCz3s4JPg9tq8wzrelgXz46qg3ItoJcd2ebq5be2dXyQ8XMLF9O8AVW6LHpE8Y0MKh/n5yeqxBTvhZzzvmeVC0/VMzM8uUEn4fWYWvZNK9pKejysACXnTg655L5hDENPHHR0fzxmk/mNU99tSz8Ui0/VDqzadMmNmzYwObNm9m8eTMbNmxg06ZNAPztb39jw4YNQDJsbsOGDUS2jiNmVpWc4PMw9aTR9Nmu/T6A2dp7u7M+fE8txtLaKa119jyK+FzF1tUfKt35fErpqquuora2lmuuuYbbbruN2tparrrqKgBGjhxJbW0tzc3NjBs3jtraWpYvX17iiM2sp6iaftE3NjbGvHnzeuS5Lp71Erc//adso+S2aqiv5YmLjn5Pz3RIkk05Jc5yjLG7veBzPb8cX7uZWS4kzY+Ixmz7ijpMrpo98ofVHSZ3eLe9txLmDS+3GAsxXC/XWQTL7bWbmRWCE3yecums1dreW84dvlpLue0tb5trjIUec96TSbecPx8zs3y5DT5PnXXWymzvLdcOX5ljxduTS4zFGHPek0m3XD8fM7PucILPU7ZOXK3adkzLt2d6sTt+ZSsldzXG9q7T3THnPZl0q2XkAFRuZ0EzKzwn+DxNGNPApw9ueM98uq2JIbMaOZ9e8D0xE1tHpeFB/fvk3MmsGKXtnky6PTVKodg8e5+ZZXIbfDdk62jXXjtxV5aNhZ5pg25veVpof8W6rlynO6XtbAvvFHMu+a5+PuXInQXNLJMTfDcUs524J9qgJ48b+Z7hYa26khiyXacQpe1qSLo9yZ0FzSyTq+i7oZjtxD3RBt1aNd2eXBNDtVRxVzp3FjSzTE7w3VDMduKeaoOeMKaBhgIkhu5OkWvdV02dBc2s+5zgu6GYJdeeLBU7MVQH16SYWSZPVWtA4SeqMTOz4vNUtdYpd2gzM6surqI3MzOrQk7wZmZmVcgJ3szMrAo5wZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRVygjczM6tCTvBmZmZVyAnezMysCjnBm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqkBO8mZlZFVJElDqGgpG0Gljew0+7C/BGDz9ntfF72H1+D7vP72Fh+H3svq68h8MiYnC2HVWV4EtB0ryIaCx1HJXM72H3+T3sPr+HheH3sfsK9R66it7MzKwKOcGbmZlVISf47rup1AFUAb+H3ef3sPv8HhaG38fuK8h76DZ4MzOzKuQSvJmZWRVygs+TpOMkLZa0RNJFpY6n0kjaQ9IjkhZJWiDpglLHVKkk1UhqkvSbUsdSqSTVS5op6Q/pd/IjpY6p0kj6Svp/+WVJd0jqV+qYKoGkX0haJenljG07S3pQ0ivpv4PyubYTfB4k1QA3AMcD+wGnSdqvtFFVnE3A1yJiX+Aw4Dy/h3m7AFhU6iAq3I+A30bEKOCD+P3sEkkNwJeBxojYH6gBTi1tVBXjFuC4NtsuAn4XESOA36WPu8wJPj+HAksiYmlEvAPMAMaXOKaKEhGvR8Rz6f23SP6gNpQ2qsojaSjwSeDmUsdSqSTtBBwB/BtARLwTEWtKG1VF2h6olbQ90B9YWeJ4KkJEPAb8pc3m8cCt6f1bgQn5XNsJPj8NwGsZj1fg5JQ3SXsBY4BnShtJRboO+BdgS6kDqWDvB1YD/542ddwsacdSB1VJIqIZ+B7wJ+B14K8R8UBpo6po74uI1yEpDAG75nMRJ/j8KMs2D0fIg6Q64G7gwohYW+p4KomkE4BVETG/1LFUuO2BDwHTImIM8DZ5Von2Vmkb8XhgOLA7sKOk00sblTnB52cFsEfG46G4OqrLJPUhSe63R8Q9pY6nAo0FTpK0jKSZ6GhJt5U2pIq0AlgREa01SDNJEr7l7ljgjxGxOiI2AvcAf1/imCrZnyUNAUj/XZXPRZzg8/MsMELScEl9STqTzC5xTBVFkkjaPBdFxA9KHU8liogpETE0IvYi+Q4+HBEuNXVRRPwv8JqkkemmY4CFJQypEv0JOExS//T/9jG4o2J3zAbOTO+fCfw6n4tsX7BwepGI2CTpfGAOSW/RX0TEghKHVWnGAp8HXpL0fLrtmxFxfwljst7rS8Dt6Q/2pcBZJY6nokTEM5JmAs+RjJBpwjPa5UTSHcBRwC6SVgCXAdcAd0o6m+TH08l5Xdsz2ZmZmVUfV9GbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3KwJJIen7GY+/LmlqD8fwqKTG9P79kuq7eb2jirFinaQn29l+i6R/zPOaHcaajtWm9TNRKsu2gyQ9la6S9qKkU/KJx6wUPA7erDj+BkyUdHVEvNHVkyVtHxGbChVMRHyiUNfqSD5xR0QpZjz7iqS1JFOqfhuYC+yfZdsy4IyIeEXS7sB8SXO8GI1VApfgzYpjE8lEH19pu0PSMEm/S0uEv5O0Z7r9Fkk/kPQI8J308TRJj0haKunIdO3oRZJuybjeNEnz0lLm5dmCkbRM0i6SzpH0fHr7Y/pcSPqHtKT6nKS70jUCkHRcukb674GJ7Vz7C+k59wIPSNoxjfPZdPGW8elxoyX9d/rcL0oakW5fl/4rSddLWijpPjIW2GiNP73fKOnR9P6hkp5Mn+fJjNnoMuM7MuM1N0kakM6euAvJEqe/jYgH2tn2PxHxCkBErCSZMnRwtvfBrNw4wZsVzw3A5yQNbLP9emB6RBwI3A78OGPfPsCxEfG19PEg4GiSHwr3Aj8ERgMHSDooPeZbEdEIHAgcKenA9gKKiJ9GxEHAISRzsP8gTZwXp8/7IWAe8FVJ/YCfAycCHwV26+C1fgQ4MyKOBr5FMm3uIcDHgGuVrM52DvCj9Pkb0+fP9ClgJHAA8P/IbS7zPwBHpIvEXAr8a5Zjvg6clz7vR4EWSRcCb5C898dJ+ni2bZkXkXQo0Bd4NYe4zErOVfRmRRIRayVNJykRtmTs+gjvlob/A/huxr67ImJzxuN7IyIkvQT8OSJeApC0ANgLeB74jKRJJP+fhwD7AS92Et6PSJLwvUpWpdsPeCJthu4LPAWMIllA5JX0OW8DJrVzvQcjonVN638gWQTn6+njfsCe6TW/pWQN+3tar5vhCOCO9PWvlPRwJ68BYCBwa1obEECfLMc8QfJD5vb0eVdI+lH6vk6NiKlp+/tDWbaRvvYhJJ/VmRHhpXmtIrgEb1Zc1wFnAx2tL545X/Tbbfb9Lf13S8b91sfbSxpOUkI9Jq0RuI8kobZL0heAYUBrdb5IEvRB6W2/iDg7S2wdyYxbwKczrrdnRCyKiF8CJ5H82Jkj6egs12nv+Tbx7t+rzNd3JfBIROxPUtPwntceEdcAXwRqgacljYp0ju6ImJr+G9m2AUjaieR9vTginu74bTArH07wZkWUlmrvJEnyrZ4kWf0N4HPA77vxFDuRJNe/SnofcHxHB0s6mOQHwekZJdGngbGSPpAe01/SPiTV38Ml7Z0ed1qOMc0BvtRaApY0Jv33/cDSiPgxyWpZbZsSHgNOlVSTlpg/lrFvGXBwev/TGdsHAs3p/S+085r3joiXIuI7JM0Po3J8HShZfOY/SZpU7sr1PLNy4ARvVnzfJ+m81erLwFmSXiRZUe+CfC8cES+QrNy1APgFSXV0R84HdgYeSTud3RwRq0mS4x1pTE8DoyJiA0mV/H1pJ7vlOYZ1JUlV+YuSXk4fA5wCvKxk9cBRwPQ25/0n8ArwEjCNpBd7q8uBH0l6HMhswvgucLWkJ0hWdszmQkkvS3qBpPbgv3J8HQCfIWk6+EJGR72DOjvJrBx4NTkzM7Mq5BK8mZlZFXKCNzMzq0JO8GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKrQ/w8c7M/FmDXnDwAAAABJRU5ErkJggg==\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": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "702.179214549006"
      ]
     },
     "execution_count": 7,
     "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": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Lagrange multiplier statistic', 4.893213374093967),\n",
       " ('p-value', 0.0865869050235217),\n",
       " ('f-value', 2.50371594625644),\n",
       " ('f p-value', 0.08794028782672986)]"
      ]
     },
     "execution_count": 8,
     "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": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('F statistic', 1.1002422436378148), ('p-value', 0.3820295068692508)]"
      ]
     },
     "execution_count": 9,
     "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": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('t value', -1.0796490077783811), ('p value', 0.283463924755859)]"
      ]
     },
     "execution_count": 10,
     "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.8.4rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
