{
 "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, 21 Feb 2020   Prob (F-statistic):           1.90e-08\n",
      "Time:                        13:56:39   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",
    "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": {
    "collapsed": true
   },
   "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": 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.713437811597181), ('Two-tail probability', 0.15618424580304824)]"
      ]
     },
     "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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfZyUdb3/8debZZEF1CXFG9YbyBtUUqMWzaOZaedgHRNSKe+SzrGfddRu1DhJeYNopWHlKW9OHUvxmJGhcSg9YWqQx8wE8Q6MBENhsdRkRWRJWD6/P65rcRhmd2dnZ3ZmZ9/Px2MfzHyvm/nMzLKf63tzfb+KCMzMzKy69Ct3AGZmZlZ8TvBmZmZVyAnezMysCjnBm5mZVSEneDMzsyrkBG9mZlaFnODNrEOSpkq6PX28l6R1kmqK/BorJH2omOds53UWSzqmnW3HSFpVpNeZJ+nTxThXpZP0KUn/V+44bFtO8NYtPfWHuZqln+FfJQ3OKPu0pHllDCuniHgxIoZERGu5YylERIyOiHnljsOsJzjBW68mqX+5YyiS/sAXunsSJar6/3WxWw/sbVX0/8lwgrcSknSCpCckNUv6naRD0vKLJc3K2vc/JH03fbyjpB9KeklSk6Sr2v6op82BD0v6jqTXgKmS9pH0oKS/SXpV0o8l1Wec+z2SFkl6Q9LPJP1U0lWdxZnj/fynpGuzyv5H0oXp4y+n8b4haamk47rwcU0HvpQZd9br/IOkxyS9nv77Dxnb5kn6mqSHgfXAO9Oyq9L3s07SLyTtlH42a9NzjMj6/Fem2xZKen87cYyQFJL6SzoiPXfbzwZJK9L9+qXf8/L0e7lT0jsyzvNJSS+k277a0Qcj6VZJN0m6V9KbwAclbSfpWkkvpq0f/ympLt1/Z0m/TL/P1yQ91HbRk9niJKkuPfcaSUuAsVmvG5L2zYrjqvTx0PQ1XkmP/6WkPdqJf19J89Pv7lVJP21nv4GSbk8/k+b0O9o13TYyPccbkn4t6Xq93W2yTddC1vs8TNIj6TlfSo8dkPU+z5P0HPBcWnZA+jqvpb/LH8/YfydJc9LflT8A+3T0/Vn5OMFbSUh6D/Aj4DPATsD3gTmStgN+AnxE0g7pvjXAx4E70sNnAJuAfYExwD8Bmf2ZhwPPA7sAXwMEfAMYDhwI7AlMTc89APg5cCvwjvS1P5ZnnNnuAD4hSemxQ9PYZkoaBZwPjI2I7YFxwIoufGQLgHnAl7I3pInxHuC7aYzfBu6RtFPGbp8EzgG2B15Iy05NyxtI/gg/AtySfg7PApdnHP8Y8O502x3AzyQN7CjgiHgkba4fAgwFfk/y+QJ8HpgAfIDke1kD3JC+n4OAm9LYhqfvKWdyzHA6yXe9PfB/wDXA/mnM+6bv8bJ034uAVcAwYFfgK0CuObkvJ/lc9iH5viZ1EkOmfiSf5d7AXkALcH07+14J3EfyGe0BfK+d/SYBO5L8/u4EfDY9LyTfyUJg5/R8XYm1FbggPfYI4Djg3Kx9JpD8vzpISVfRr9PX3AU4DbhR0uh03xuADcDuwL+mP1aJIsI//in4hySJfShH+U3AlVllS4EPpI//DzgrffyPwPL08a7A34G6jONOA36TPv4U8GInMU0AFqWPjwaaAGVs/z/gqnzizCoX8CJwdPr8/wEPpo/3BV4GPgTUFvIZAu8CXidJTJ8G5qXbPwn8IeuYR4BPpY/nAdOyts8Dvprx/FvA/2Y8/yjwRAcxrQEOTR9PBW5PH48gSZb9c3zf9wD90ufPAsdlbN8d2EjSFXEZMDNj22DgrVy/R+n2W4Hbsr6HN4F9MsqOAP6cPp4G/A+wb0e/ryQXicdnbDsHWJXxPDLPkcZxVTsxvhtYk/X5fzp9fBvwA2CPTn4P/hX4HXBIVvleJBe8gzPK7sj4To7JjLuj/5fpti8CP896n8dmPP8E8FDWMd8nuSCqSb/HAzK2fR34v678zvunZ35cg7dS2Ru4KG0WbJbUTFIzGZ5uv4MkcUNSO7sj47ha4KWM475PUpNoszLzhSTtImmmkubxtcDtJLUV0tdrivQvUY7jO4tzi/QcM7Pi/nG6bRnJH86pwMtpPNucoyMR8QzwS+DirE3DebtW3uYFklprrvfU5q8Zj1tyPB/S9kTSRZKeTZuRm0lqkjuTB0mfIUkyp0fE5rR4b+DnGZ/psyQ1yV3T97Ml3oh4E/hbJy+T+f6GAYOAhRnn/1VaDkl3xzLgPknPS8r+PNtsFQfbfsbtkjRI0vfTboa1wG+BeuUeH/DvJBclf1Ayir+9Gu9/A3NJWoRWS/qmpNo0zjXp51RIrPunXQh/SWP9Ott+t9n/Jw7P+j9xBrAbyWfcnwI/N+tZTvBWKiuBr0VEfcbPoIhoa8L9GXBM2m/5Md5O8CtJavA7Zxy3Q0SMzjh3dnPrN9KyQyJiB+BMkj+oAC8BDW3N6qk9uxBntp8Ap0jam6RJ864tQUXcERFHkfyBDJJm5K66nKRlIDN5r07PmWkvkpaJLS9fwGsBoKS//csk3SRDI6KepCVBHR749rFXAuMj4vWMTSuBD2d9rgMjoonkO9kz4xyDSJqkO5L5/l4luUAZnXHuHSPpKiAi3oiIiyLinSQtFRcq93iIreIg+UwzrSe5kGizW8bji4BRwOHp79zRbW9nm8Aj/hIR/y8ihpN0Bd2Y2befsd/GiLgiIg4C/gE4ATgrjXOoMu6yyIr1zcw404uMYRnbbwL+COyXxvqVHHFmXwDPz/ruhkTEvwGvkLQmdPS5WYVwgrdiqE0HCLX99Af+C/ispMOVGCzpnyVtDxARr5A0Y95C0rT6bFr+Ekl/5bck7aBksNY+kj7QwetvD6wDmiU1AJMztj1CUnM8X8nAsPHAYRnbO4wzW0QsIvkjdzMwNyKaASSNknRs2ne/gSQBdflWsrQl4Kckfdht7gX2l3R6+h4+ARxEUtsvhu1J/mi/AvSXdBmwQ2cHSdozjfWsiPhT1ub/BL6WXgghaVj62QPMAk6QdFQ6RmIaXfhblLYS/BfwHUm7pOdvkDQufXyCkoFtAtaSfA+5vos7gSlKBsztAXwua/sTwOmSaiQdTzKeoM32JN9xczpG4nLaIWmi3h6At4YkmW4Tj6QPSjo4TdBrSZrCWyPiBZIxGldIGiDpKJILlzZ/Agamv7e1wCVA5hiS7dPzrZN0APBv7cWa+iXJ79snJdWmP2MlHRjJ7ZF3kwxuHZSOp+jKeADrQU7wVgz3kvyxa/uZGhELSGqi15P8UVtG0n+e6Q6Svuc7ssrPAgYAS9JjZ5H04bbnCuA9JLXOe0j+AAEQEW8BJwFnA80ktftfkrQSkGec2X6SI+7tgKtJapd/IelS+AqApDMkLe7knJmmkfRLt72Hv5HU5i4iacr+d+CEiHi1C+fsyFzgf0kSxQskFyi5mvyzHUdSq52lt0fSt73P/wDmkDSTv0EyAO/w9P0sBs4j+fxeIvncuzrBzJdJvqvfp83O95PUqAH2S5+vI7nAuzFy3/t+Bcn7/TPJReV/Z23/AkkibWuinp2x7TqgjuT7/j1JF0F7xgKPSlpH8pl8ISL+nGO/3Uh+19eSdGnMJ+lugqQ76HDgNZKLidvaDkpbTs4luehsIqnRZ36eX0qPf4PkwijnKP6M871BMnj0VJLWo7+QtEa1XTScT9K98xeScQm3dHQ+Kx9t3TVpVv0kPQr8Z0T4D5P1SpKmkgwAPLPcsVjlcg3eqp6kD0jaLW3engQcQsc1LjOzXs+zFllfMIqkv3UIsBw4Je3rNzOrWm6iNzMzq0JuojczM6tCTvBmZmZVqGr64HfeeecYMWJEucMwMzPrMQsXLnw1Iobl2lY1CX7EiBEsWLCg3GGYmZn1GEntThXsJnozM7Mq5ARvZmZWhZzgq8Rzzz3HwIEDOfNMT2xlZmZO8FXjvPPOY+zYseUOw8zMKoQTfBWYOXMm9fX1HHdcrhUxzcysL3KC7+XWrl3LZZddxre+9a1yh2JmZhXECb6Xu/TSSzn77LPZc889yx2KmZlVkKq5D74veuKJJ7j//vtZtGhRuUMxM7MK4wTfi82bN48VK1aw1157AbBu3TpaW1tZsmQJjz/+eJmjMzOzcqqa1eQaGxujr81kt379etauXbvl+bXXXsuKFSu46aabGDYs58yFZmZWRSQtjIjGXNtcg+/FBg0axKBBg7Y8HzJkCAMHDnRyNzMzJ/hqMnXq1HKHYGZmFcKj6M3MzKqQE7yZmVkVKmmCl3S8pKWSlkm6OMf2oyU9LmmTpFNybN9BUpOk60sZp5mZWbUpWYKXVAPcAHwYOAg4TdJBWbu9CHwKuKOd01wJzC9VjGZmZtWqlIPsDgOWRcTzAJJmAuOBJW07RMSKdNvm7IMlvRfYFfgVkPMWAIPZi5qYPncpq5tbGF5fx+Rxo5gwpqHcYZmZWZmVsom+AViZ8XxVWtYpSf2AbwGTO9nvHEkLJC145ZVXCg60t5q9qIkpdz9NU3MLATQ1tzDl7qeZvaip3KGZmVmZlTLBK0dZvrPqnAvcGxErO9opIn4QEY0R0dgX7/2ePncpLRtbtypr2djK9LlLyxSRmZlVilI20a8CMldA2QNYneexRwDvl3QuMAQYIGldRGwzUK8vW93c0qVyMzPrO0qZ4B8D9pM0EmgCTgVOz+fAiDij7bGkTwGNTu7bGl5fR1OOZD68vq4M0ZiZWSUpWRN9RGwCzgfmAs8Cd0bEYknTJJ0IIGmspFXAROD7khaXKp5qNHncKOpqa7Yqq6utYfK4UWWKyMzMKoUXm+nlPIrezKzv8mIzVWzCmAYndDMz24anqjUzM6tCTvBmZmZVyAnezMysCjnBm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqkBN8noYMGbLVT01NDZ/73Oe2bH/ggQc44IADGDRoEB/84Ad54YUXyhitmZn1dU7weVq3bt2Wn7/+9a/U1dUxceJEAF599VVOOukkrrzySl577TUaGxv5xCc+UeaIzcysL3OCL8CsWbPYZZddeP/73w/A3XffzejRo5k4cSIDBw5k6tSpPPnkk/zxj38sc6RmZtZXOcEXYMaMGZx11llIyZL3ixcv5tBDD92yffDgweyzzz4sXuy1c8zMrDyc4LvoxRdfZP78+UyaNGlL2bp169hxxx232m/HHXfkjTfe6OnwzMzMACf4Lrvttts46qijGDly5JayIUOGsHbt2q32W7t2Ldtvv31Ph2dmZgY4wXfZbbfdtlXtHWD06NE8+eSTW56/+eabLF++nNGjR/d0eGZmZoATfJf87ne/o6mpacvo+TYf+9jHeOaZZ7jrrrvYsGED06ZN45BDDuGAAw4oU6RmZtbXOcF3wYwZMzjppJO2aXofNmwYd911F1/96lcZOnQojz76KDNnzixTlGZmZqCIKHcMRdHY2BgLFiwodxhmZmY9RtLCiGjMtc01eDMzsyrkBG9mZlaFnODNzMyqkBO8mZlZFepf7gB6s9mLmpg+dymrm1sYXl/H5HGjmDCmodxhmZmZOcEXavaiJqbc/TQtG1sBaGpuYcrdTwM4yZuZWdm5ib5A0+cu3ZLc27RsbGX63KVlisjMzOxtTvAFWt3c0qVyMzOznuQEX6Dh9XVdKjczM+tJTvAFmjxuFHW1NVuV1dXWMHncqDJFZGZm9jYPsitQ20A6j6I3M7NK5ATfDRPGNDihm5lZRXITvZmZWRVygjczM6tCTvBmZmZVyAnezMysCjnBm5mZVaGSJnhJx0taKmmZpItzbD9a0uOSNkk6JaP83ZIekbRY0lOSPlHKOM3MzKpNyRK8pBrgBuDDwEHAaZIOytrtReBTwB1Z5euBsyJiNHA8cJ2k+lLFamZmVm1KeR/8YcCyiHgeQNJMYDywpG2HiFiRbtuceWBE/Cnj8WpJLwPDgOYSxmtmZlY1StlE3wCszHi+Ki3rEkmHAQOA5UWKy8zMrOqVMsErR1l06QTS7sB/A/8SEZtzbD9H0gJJC1555ZUCwzQzM6s+pUzwq4A9M57vAazO92BJOwD3AJdExO9z7RMRP4iIxohoHDZsWLeCNTMzqyalTPCPAftJGilpAHAqMCefA9P9fw7cFhE/K2GMZmZmValkCT4iNgHnA3OBZ4E7I2KxpGmSTgSQNFbSKmAi8H1Ji9PDPw4cDXxK0hPpz7tLFauZmVm1UUSXusUrVmNjYyxYsKDcYZiZmfUYSQsjojHXNs9kZ2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7wZmZmVcgJ3szMrAo5wZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRVygjczM6tCTvBmZmZVyAnezMysCjnBm5mZVSEneDMzsyrkBG9mZlaFnODzdMwxxzBw4ECGDBnCkCFDGDVq1JZtr7zyCqeffjr19fUMHTqUM844o4yRmpmZQf9yB9CbXH/99Xz605/epvykk05i7NixvPDCCwwaNIhnnnmmDNGZmZm9zQm+m+677z5WrlzJvHnzqKmpAWDMmDFljsrMzPo6N9F3wZQpU9h555058sgjmTdvHgC///3vGTVqFJMmTWKnnXZi7NixzJ8/v7yBmplZn+cEn6drrrmG559/nqamJs455xw++tGPsnz5clatWsV9993HBz/4Qf7yl79w0UUXMX78eF599dVyh2xmZn2YE3yeDj/8cLbffnu22247Jk2axJFHHsm9995LXV0dI0aM4Oyzz6a2tpZTTz2VPffck4cffrjcIZuZWR/mBF8gSUQEhxxyCJLKHY6ZmdlWnODz0NzczNy5c9mwYQObNm3ixz/+Mb/97W8ZN24cH/vYx1izZg0zZsygtbWVWbNm0dTUxJFHHlnusM3MrA/zKPo8bNy4kUsuuYQ//vGP1NTUcMABBzB79uwt98LPmTOHc889l/POO48DDjiA//mf/2HnnXcuc9RmZtaXKSLKHUNRNDY2xoIFC8odhpmZWY+RtDAiGnNtcxO9mZlZFXKCNzMzq0JO8GZmZlXIg+wKNHtRE9PnLmV1cwvD6+uYPG4UE8Y0lDssMzMzwAm+ILMXNTHl7qdp2dgKQFNzC1PufhrASd7MzCqCm+gLMH3u0i3JvU3Lxlamz11apojMzMy25gRfgNXNLV0qNzMz62lO8AUYXl/XpXIzM7Oe5gRfgMnjRlFXW7NVWV1tDZPHjSpTRGZmZlsraYKXdLykpZKWSbo4x/ajJT0uaZOkU7K2TZL0XPozqZRxdtWEMQ1846SDaaivQ0BDfR3fOOlgD7AzM7OKUbJR9JJqgBuAfwRWAY9JmhMRSzJ2exH4FPClrGPfAVwONAIBLEyPXVOqeLtqwpgGJ3QzM6tYpazBHwYsi4jnI+ItYCYwPnOHiFgREU8Bm7OOHQf8OiJeS5P6r4HjSxirmZlZVSllgm8AVmY8X5WWFe1YSedIWiBpwSuvvFJwoGZmZtWmlAleOcryXbour2Mj4gcR0RgRjcOGDetScGZmZtWslAl+FbBnxvM9gNU9cKyZmVmfV8oE/xiwn6SRkgYApwJz8jx2LvBPkoZKGgr8U1pmZmZmeShZgo+ITcD5JIn5WeDOiFgsaZqkEwEkjZW0CpgIfF/S4vTY14ArSS4SHgOmpWVmZmaWB0Xk2y1e2RobG2PBggXlDsPMzKzHSFoYEY25tnkmOzMzsyrkBG9mZlaFnODNzMyqkBN8ETz33HMMHDiQM888E4B58+bRr18/hgwZsuVnxowZZY7SzMz6kpLNRd+XnHfeeYwdO3arsuHDh7Nq1aoyRWRmZn2da/DdNHPmTOrr6znuuOPKHYqZmdkWTvDdsHbtWi677DK+9a1vbbPt5ZdfZtddd2XkyJFccMEFvPnmm2WI0MzM+ion+G649NJLOfvss9lzzz23Kj/ggAN44okneOmll3jwwQdZuHAhF154YZmiNDOzvsgJvkBPPPEE999/PxdccME223bbbTcOOugg+vXrx8iRI/nmN7/JrFmzyhClmZn1VR5kV6B58+axYsUK9tprLwDWrVtHa2srS5Ys4fHHH99qX0lUy4yBZmbWO3iq2gKtX7+etWvXbnl+7bXXsmLFCm666SYWL17MO9/5Tvbcc09WrVrFWWedxYgRI7jlllt6LD4zM6t+HU1V6xp8gQYNGsSgQYO2PB8yZAgDBw5k2LBhPP7445xxxhmsWbOGnXbaiQkTJvD1r3+9jNGamVlf4xq8mZlZL+XFZszMzPoYJ3gzM7Mq5ARvZmZWhTzIrohmL2pi+tylrG5uYXh9HZPHjWLCmIZyh2VmZn2QE3yRzF7UxJS7n6ZlYysATc0tTLn7aQAneTMz63Fuoi+S6XOXbknubVo2tjJ97tIyRWRmZn2ZE3yRrG5u6VK5mZlZKTnBF8nw+roulZuZmZWSE3yRTB43irramq3K6mprmDxuVJkiMjOzvsyD7IqkbSCdR9GbmVklcIIvogljGpzQzcysIriJ3szMrAo5wZuZmVUhJ3gzM7Mq5ARvZmZWhfJK8JJ2lfRDSf+bPj9I0tmlDc3MzMwKlW8N/lZgLjA8ff4n4IulCMjMzMy6L98Ev3NE3AlsBoiITUBrx4eYmZlZueSb4N+UtBMQAJLeB7xesqjMzMysW/Kd6OZCYA6wj6SHgWHAKSWLyszMzLolrwQfEY9L+gAwChCwNCI2ljQyMzMzK1heCV7SSVlF+0t6HXg6Il4uflhmZmbWHfn2wZ8N3Ayckf78F0mz/cOSPtneQZKOl7RU0jJJF+fYvp2kn6bbH5U0Ii2vlTRD0tOSnpU0pYvvy8zMrE/LN8FvBg6MiJMj4mTgIODvwOHAl3MdIKkGuAH4cLr/aZIOytrtbGBNROwLfAe4Ji2fCGwXEQcD7wU+05b8zczMrHP5JvgREfHXjOcvA/tHxGtAe33xhwHLIuL5iHgLmAmMz9pnPDAjfTwLOE6SSEbrD5bUH6gD3gLW5hmrmZlZn5fvKPqHJP0S+Fn6/GTgt5IGA83tHNMArMx4voqkxp9zn4jYlPbr70SS7McDLwGDgAvSiwkzMzPLQ74J/jySpH4kySj624C7IiKAD7ZzjHKURZ77HEYykc5wYCjJBcb9EfH8VgdL5wDnAOy11175vRMzM7M+IK8m+kjMiogLIuKL6ePsZJ1tFbBnxvM9gNXt7ZM2x+8IvAacDvwqIjamo/QfBhpzxPWDiGiMiMZhw4bl81Yq1plnnsnuu+/ODjvswP7778/NN99c7pDMzKwXy3exmfdJekzSOklvSWqV1Fmf+GPAfpJGShoAnEoyWU6mOcCk9PEpwIPphcOLwLFKDAbeB/wx3zfVG02ZMoUVK1awdu1a5syZwyWXXMLChQvLHZaZmfVS+Q6yux44DXiOZNDbp4HvdXRAOl/9+SSL1DwL3BkRiyVNk3RiutsPgZ0kLSO57a7tVrobgCHAMyQXCrdExFN5v6teaPTo0Wy33XYASEISy5cvL3NUZmbWW+XbB09ELJNUExGtwC2SfpfHMfcC92aVXZbxeAPJLXHZx63LVV7tzj33XG699VZaWloYM2YMH/nIR8odkpmZ9VL51uDXp83sT0j6pqQLgMEljKtPuvHGG3njjTd46KGHOOmkk7bU6M3MzLoq3wT/yXTf84E3SQbGnVyqoPqympoajjrqKFatWsVNN91U7nDMzCyHmTNncuCBBzJ48GD22WcfHnroIZYsWUJjYyNDhw5l6NChfOhDH2LJkiVli7HTJvp0RrqvRcSZwAbgipJHZWzatMl98GZmFejXv/41X/7yl/npT3/KYYcdxksvvQTA4MGDmTVrFnvvvTebN2/mhhtu4NRTT+Wpp8ozhKzTGnza5z4sbaK3Enj55ZeZOXMm69ato7W1lblz5/KTn/yEY489ttyhmZlZlssvv5zLLruM973vffTr14+GhgYaGhqor69nxIgRSCIiqKmpYdmyZWWLM99BditIFpaZQ9JED0BEfLsUQfU1krjpppv47Gc/y+bNm9l777257rrrGD8+e2ZfMzMrp9bWVhYsWMCJJ57Ivvvuy4YNG5gwYQLTp0+nrq4OgPr6etatW8fmzZuZNm1a2WLNN8GvTn/6AduXLpy+adiwYcyfP7/cYZiZWSf++te/snHjRmbNmsVDDz1EbW0t48eP56qrruJrX/saAM3Nzbz55pvMmDGDvffeu2yxqvMJ6TJ2lgZHxJud79nzGhsbY8GCBeUOw8zMqtiaNWt4xzvewa233sqkSck8bXfddRdXXXUVixYt2mrfzZs3M2zYMJ599ll22WWXksQjaWFEbDPTK+Q/k90RkpaQTFiDpEMl3VjEGM3MzCre0KFD2WOPPUgWPu3Y5s2bWb9+PU1NTT0Q2bbyvU3uOmAc8DeAiHgSOLpUQZmZmVWqf/mXf+F73/seL7/8MmvWrOG6667jhBNO4Ne//jWLFi2itbWVtWvXcuGFFzJ06FAOPPDAssSZb4InIlZmFbUWOZY+bfaiJo68+kFGXnwPR179ILMXleeKz8zMOnbppZcyduxY9t9/fw488EDGjBnDV7/6VZqbmznttNPYcccd2WeffVi2bBm/+tWvGDhwYFnizKsPXtIs4Nskc9K/D/g80BgRp5Y2vPz15j742YuamHL307RsfPuaqa62hm+cdDATxjSUMTIzM6tk3e6DBz5LsiZ8A8kSr+9On1sRTJ+7dKvkDtCysZXpc5eWKSIzM+vt8r1NThFxRkkj6cNWN7d0qdzMzKwz+dbgfyfpPklnS6ovaUR90PD6ui6Vm5mZdSavBB8R+wGXAKOBxyX9UtKZJY2sD5k8bhR1tTVbldXV1jB53KgyRWRmZr1dV9aD/wPwB0lfJxlwNwO4vVSB9SVtA+mmz13K6uYWhtfXMXncKA+wMzPrBWYvaqrIv995JXhJOwAfA04F9gF+DhxWwrj6nAljGiriF8LMzPKXfRdUU3MLU+5+GqDsf9Pz7YN/kmTk/LSI2D8ivhwRC0sYl5mZWcWr5Lug8m2if2dEhKTBJY3GzMysF6nku6DyrcG/z3PRm5mZba2S74LyXPRmZmYFquS7oLoyin5l1uo5novezMz6tEq+CyrfBL9S0j8AIWkAyVz0z5YuLOtIpU/fmGoAAB4kSURBVN6SYWbWF1XqXVD5JvjPAv/B23PR3wecW6qgrH2VfEuGmZlVjnxnsns1Is6IiF0jYpeIOBM4q8SxWQ6VfEuGmZlVjrzXg8/hwqJFYXlr79aLpgq4JcPMzCpHdxK8Ot/Fiq29Wy9E0nxvZmYG3UvwUbQoLG+Tx43KeWUV4GZ6MzPbosMEL+kNSWtz/LwBDO+hGC3DhDEN7V5ZVcLMSWZmVhk6HEUfEdv3VCCWv4b6upx97pUwc5KZmVWG7jTRW5lU8sxJZmZWGfKeyc4qRyXPnGRmZpXBCb6XqtSZk8zMrDK4id7MzKwKuQZvPcrz6JuZ9QwneOsxnkffzKznlLSJXtLxkpZKWibp4hzbt5P003T7o5JGZGw7RNIjkhZLelrSwFLGaqXnefTNzHpOyRK8pBrgBuDDwEHAaZIOytrtbGBNROwLfAe4Jj22P3A78NmIGA0cA2wsVazWM9qbiMcT9JiZFV8pm+gPA5ZFxPMAkmYC44ElGfuMB6amj2cB10sS8E/AUxHxJEBE/K2EcVoPGV6BE/R4TICZVatSNtE3ACsznq9Ky3LuExGbgNeBnYD9gZA0V9Ljkv491wtIOkfSAkkLXnnllaK/ASuuSpugp21MQFNzC8HbYwK8aI+ZVYNS1uDbWxMln336A0cBY4H1wAOSFkbEA1vtGPED4AcAjY2NFbH4jWuE7au0CXo6GhPg78zMertSJvhVwJ4Zz/cAVrezz6q0331H4LW0fH5EvAog6V7gPcADVDCPEu9cJU3Q4zEBZlbNStlE/xiwn6SRkgYApwJzsvaZA0xKH58CPBgRAcwFDpE0KE38H2DrvvuK5FHivUt7ff9etMfMqkHJEnzap34+SbJ+FrgzIhZLmibpxHS3HwI7SVoGXAhcnB67Bvg2yUXCE8DjEXFPqWItFtcIe5dKGxNgZlZMJZ3oJiLuBe7NKrss4/EGYGI7x95Ocqtcr1GJo8StfZU2JsDMrJg8k10RTR43aqs+eHCNsNJV0pgAM7NicoIvItcIzcysUjjBF5lrhGZmVgm8XKyZmVkVcoI3MzOrQk7wZmZmVcgJ3szMrAp5kF0ZeL56MzMrNSf4Hub56s3MrCe4ib6Heb56MzPrCU7wPczz1ZuZWU9wgu9hXsHMzMx6ghN8D5q9qIn1b23aptzz1ZuZWbF5kF0PyR5c16a+rpapJ472ADszMysqJ/gi6uj2t1yD6wAGb9d/yz6+fc7MzIrFCb5IOrv9rbPBdb59zszMisl98EXS2e1vnQ2u8+1zZmZWTE7wRdJZDX3yuFHU1dZstS1zcJ1vnzMzs2Jygi+SzmroE8Y08I2TDqahvg4BDfV1fOOkg7c0v/v2OTMzKyb3wRfJ5HGjthkln33724QxDe32p+dzvJmZWb5cgy+SCWMaOPm9DdRIANRInPze9hN6ruM7quGbmZl1hWvwRTJ7URN3LWyiNQKA1gjuWthE497v6FKS720J3bf2mZlVJif4bmpLcE05BsO1jYKv9IRXaJL2rX1mZpXLTfTd0JbgciX3NpU+Cj7zPQRvJ+nZi5o6Pda39pmZVS7X4Lth6pzFOWeny5TvKPhyNXV3lKQ7e33f2mdmVrlcgy/Q7EVNNLds7HCffEfBd6cW3V3dSdK+tc/MrHI5wReos2boroyCL2dTd3eSdGeT95iZWfm4ib5AHdVwr/vEu7vUvF7Opu7u3H+fuZCOR9GbmVUWJ/gCDa+vyzm4buig2i4nuPbO1RNN3d1N0r3x1j4zs77ACb5A7dV8L//o6KKdq6eaup2kzcyqjxN8gYrZPN2TTd2emMbMrG9QpDOv9XaNjY2xYMGCcodR0bInpoGkpcBT4pqZ9U6SFkZEY65trsEXWSXXkLtzz7uZmfUuTvBF1BNTt3bnAsIT05iZ9R1O8N2QnWzXv7WppDXk7l5AlHO0vpmZ9aySTnQj6XhJSyUtk3Rxju3bSfppuv1RSSOytu8laZ2kL5UyzkLMXtTE5FlPbjX73Jr1uWe2K1YNubsT4nhiGjOzvqNkCV5SDXAD8GHgIOA0SQdl7XY2sCYi9gW+A1yTtf07wP+WKsbuuOIXi9nYmt8AxWLVkLvbxO41583M+o5SNtEfBiyLiOcBJM0ExgNLMvYZD0xNH88CrpekiAhJE4DngTdLGGPB2qutZytmDbl+UG3O1+3KBYTveTcz6xtKmeAbgJUZz1cBh7e3T0RskvQ6sJOkFuDLwD8CFdc835mG+rqij6KfvaiJdRs2bVNeWyM3seepku9wMDMrtlImeOUoy27Tbm+fK4DvRMQ6Kdcu6cHSOcA5AHvttVeBYRamvq4252py9XW1PHzxsUV/velzl7Jx87ZdAoMH9HeSykNP3OFgZlZJSjnIbhWwZ8bzPYDV7e0jqT+wI/AaSU3/m5JWAF8EviLp/OwXiIgfRERjRDQOGzas+O+gA1NPHE1tv60vPmr7iakndn2q2ny018/+eidL1lqinCv2mZmVQylr8I8B+0kaCTQBpwKnZ+0zB5gEPAKcAjwYydR672/bQdJUYF1EXF/CWLusK9PLFqNp2Le4dY/nADCzvqZkNfiI2AScD8wFngXujIjFkqZJOjHd7Yckfe7LgAuBbW6lq1T5Ju22puHM2+mm3P00sxc1den1fItb93Rn3Xszs97Ic9EXoCtzuh959YM5a94N9XVd7qv3ILHCeR5+M6tGnou+yLoyp3sxm4Z9i1vhenLFPjOzSuAEX4CuJG33nVcOXyCZWV9S0qlqq1VX+nP7Ut/57EVNHHn1g4y8+B6OvPrBLo8zMDOz4nGCL0BXknZfmR62WIMJzcysONxEX4Cu9uf2habh3rjWvActmlk1c4IvUCFJu5oTSm+7z9wz25lZtXMTfQ/Jtbzs5FlPVk0Tdm+7z9wz25lZtXOC7yG5lpfd2Bpc8YvFZYqouHrbYMLe1uJgZtZVbqLvIe0tL5tZ3pub8Hvbfea+fdHMqp0TfIWohj7h3jSYcPK4UTlntqvUFgczs65yE30Pqa+r7bDcfcI9q6/cvmhmfZdr8D1k6omjmfyzJ7da0z1zeVn3Cfe83tTiYGbWVa7BF6irs7ZNGNPA9ImHblVjnD7x0C0Jpr2+3/pBtZ4dzqybhgwZstVPTU0Nn/vc57Zsv/nmm9l3330ZMmQIxx9/PKtXry5jtGbF4dXkClCKlclynbO2RhBsVev3Cmhm3fPmm2+y6667cu+993L00Uczf/58Jk6cyG9+8xv2228/vvCFL7BkyRLmz59f7lDNOtXRanKuwRegFP3lufqEBw/ov1VyL8brmPV1s2bNYpddduH9738/AL/4xS+YOHEio0ePZsCAAVx66aX89re/Zfny5WWO1Kx73AdfgPb6xXPddtUV2X3CIy++J6/X782315n1tBkzZnDWWWchCYCIILMls+3xM888wz777FOWGM2KwTX4ArTXXy4oah95PrPDlWORF68aZ73Viy++yPz585k0adKWso985CPceeedPPXUU7S0tDBt2jQksX79+jJGatZ9TvAFaO9e6YCiNp/nMztcT99e5wsK681uu+02jjrqKEaOHLml7LjjjuOKK67g5JNPZu+992bEiBFsv/327LHHHmWM1Kz7nOCLrLPb2rqSrPK5V7snb6+bvaiJi+58suovKKx63XbbbVvV3tucd955PPfcc7z88sucfPLJbNq0iXe9611liNCseNwHX4Cpc9qfP76fxMiL79nSFw5vT99aP6iWdRs2bRk4l89sdZ3dq91TU662JdrWdu66KNX9+r1xGVqrTL/73e9oampi4sSJW5Vv2LCBZcuWMXr0aFauXMk555zDF77wBYYOHVqmSM2KwzX4AjS35J5XHqA14u3V4n725FYryK1Zv7Hoo+J7apGXXIk2U6nmcG9v4KInALKumjFjBieddBLbb7/9VuUbNmzg9NNPZ8iQIRx22GEcccQRXHnllWWK0qx4XIMvoexk3p6m5haOvPrBgkbBF7rIS1dH3neUUEs1h/vsRU2IZGxDNi8KY131/e9/P2d5fX09Tz31VA9HY1Z6TvAFGDqott3V4QrVVlMtZJGZrk65WsjCNu11BdRIJZt4Z/rcpTmTu2h/oKOZmSXcRF+Ayz86uqTnL/VkNoWMvG+vK+BbHz+0ZH3h7bUaBL1nhb1iOvPMM9l9993ZYYcd2H///bn55pu32eeKK65AEvfff38ZIjSzSuIEX4AJYxoYVNv5R1fbT8l0swUoZR9zISPvy7H6WnvN8A19tHl+ypQprFixgrVr1zJnzhwuueQSFi5cuGX78uXLmTVrFrvvvnsZozSzSuEm+gIN6F/D+o2btylX2mmcaxR9e83cuZSyj7nQkfc9vfqa12zf2ujRb7ccSUISy5cv573vfS8A559/Ptdccw3nnntuuULsVTwDpFU7J/gCvd7eSPqAP1/9z1sVZf7ROPLqB/NK8qVMYr0lcRY6gLCanXvuudx66620tLQwZswYPvKRjwDws5/9jAEDBmx5bh0rZByKWW/jBF+gQmvBuZJrtqGDaovyR6a9GkpvSpxes31rN954I9/73vd45JFHmDdvHttttx3r1q3jK1/5Cvfdd1+5w+s1PL+C9QVO8AUqtBacmVybmlu2uQ2srramKIP4OquhOHH2XjU1NRx11FHcfvvt3HTTTbzwwgt88pOf3Gr6VetYT84AaVYuHmRXoO4MOpswpoGHLz6WFVf/M9/5xLu7dI58p7rt6Tnqredt2rSJ5cuX88ADD/Dd736X3Xbbjd12242VK1fy8Y9/nGuuuabcIVasfBZyMuvtXIPvhmLUgrtyjq70G/a2GooHPHXs5Zdf5sEHH+SEE06grq6O+++/n5/85CfccccdXHbZZWzc+PaYkLFjx/Ltb3+bD3/4w2WMuLL1lnEoZt3hBN+LdKXfsBRz1JcqCXvAU+ckcdNNN/HZz36WzZs3s/fee3Pdddcxfvz4bfatqalh6NChDBkypAyR9g69aRyKWaGc4HuRrtTKi11DKWUS9oCnzg0bNoz58+fnte+KFStKG0yV8DgUq3bugy+RUqxh3pV+w2JPTFPKPv3e1p1gZtYbuAZfArlqu1/86RNccOcTnHH4Xlw14eCCztvVWnkxayilTMI9teStmVlf4hp8CbS3tGoE3P77F7lk9tNdPmdb/3fLxlZqlEx/2xPTxbYp5ajjnlry1sysLylpgpd0vKSlkpZJujjH9u0k/TTd/qikEWn5P0paKOnp9N9jSxlnsXVWq/3Joyu7dL62FoG2Wm5rxJYE2FN9iKVMwuWY597MrNqVrIleUg1wA/CPwCrgMUlzImJJxm5nA2siYl9JpwLXAJ8AXgU+GhGrJb0LmAv0mr/2nc053xrtrxOfa6R6JQxCK/WoYw946jrfWmhmHSllH/xhwLKIeB5A0kxgPJCZ4McDU9PHs4DrJSkiFmXssxgYKGm7iPh7CeMtms6mo21rYs/W3kj19s7T04PQnIQrh28tNLPOlDLBNwCZbdGrgMPb2yciNkl6HdiJpAbf5mRgUa7kLukc4ByAvfbaq3iRd0FHtaiv3P1UzhXnTjt8z5znaq+mXiPlrPV7EFrfVQmtOmZW2UrZB5+rmpqdpTrcR9Jokmb7z+R6gYj4QUQ0RkTjsGHDCg60UJl948HbtajZi5qYMKaBJVd+mDPft9eWGnuNxJnva38UfXs18rY+90wehNa3+dZCM+tMKWvwq4DMquoewOp29lklqT+wI/AagKQ9gJ8DZ0XE8hLGWbB8alFXTTg479vi2uu7r6+r3XJuSFabu/yjo7epqblPtu9o+11Zu/AXvPnMA7z1ygoGH/gBDj3jKwC89dZbnH766SxYsIAXXniB3/zmNxxzzDHlDdrMelQpa/CPAftJGilpAHAqMCdrnznApPTxKcCDERGS6oF7gCkR8XAJY+yWYteico1Ur+0n3nxrE80Z689vyNHs31FrglWftt+V/kN2YscjPsGQg/+Rmn7aqlWnbcW53XbbrYyRmlm5lKwGn/apn08yAr4G+FFELJY0DVgQEXOAHwL/LWkZSc391PTw84F9gUslXZqW/VNEvFyqeLtq9qImtlnrNVVo33iukerr39rEmvUbt9qvZWMrU+cs3qp2Xqw+WbcC9A5bflcGD2B1cwvbNa9gvyFvbSkfMGAAX/ziF4Fkbnoz63tKOpNdRNwL3JtVdlnG4w3AxBzHXQVcVcrYumP2oiYu+tmT5LrbrTarFtVV2SPVR158T879mls2bunrh+K0JlTbyOxqv1jJ/F255JJHWLVqVZkjMrNK4pnsCnDFLxbTujn3vexDBvYvahLpqDUgcx74Ysw0V01ryLvLwsz6Oif4AmQ3mWdq7mBbITpqDcisnRdjprlqGpldTRcrZmaFcIIvsmLfmz5hTANDB9V2+lrFmO61lPPN97RqulgxMyuEE3wB2m5by6WpuaVoy8O2ufyjo/OqnU8Y08DDFx/Ln6/+Zx6++NgudxVU06Iv1XSx0pFNmzaxYcMGWltbaW1tZcOGDWzatAmAv//972zYsAFIbpvbsGED0cE0yWZWXZzgCzD1xNHU9ss93Szk7u/tzvrwPbkYy8Dat38l6utqe+2iL129WOnO91NOV111FXV1dVx99dXcfvvt1NXVcdVVyfjUUaNGUVdXR1NTE+PGjaOuro4XXnihzBGbWU9RtVzRNzY2xoIFC3rs9S6Z/TQ//v2Lue6S26Khvo6HLz52m9HpkCSbSkqelRhjd0fB53t8Jb53M7N8SFoYEY25tpX0Nrlq9ps/vtJhcoe3+3t7w7zhlRZjMW7Zy3dxnEp772ZmxeAEX6B8Bmu19fdW8oCvtlpue8vb5htjse8578mkW8nfj5lZodwHX6DOBmtl9vdW6oCvzHvF25NPjKW457wnk26lfj9mZt3hBF+gyeNGUVuTe6Bd9iC4Qkan98Sgr1y15K7E2NF5unvPeU8m3Wq6e6C3DhY0s+JzE313ZHXC1/YT0yceuk0Tcq455jtqwu6pKWM7qg23t2JdV87Tndr25HGjcg58K0XS7er3U6mqbaphM+seJ/gCTZ+7lI1Z09Vu3Bzt9hHnO+Cr7dw90f/c3vK0kHvFuq6epzu17Z5Oul35fiqVBwuaWSYn+AKVso+4p/qfc9WS23QlMZSqtl0NSbcnebCgmWVyH3yBStlH3FP9z20T6LQn38TQkxPxWPs8WNDMMjnBF6iUA7N6ctDXhDENNBQhMXR3mlzrvmoaLGhm3ecEX6BS1lp7ukbsxFAd3JJiZpk8Va0BxZ+oxszMSs9T1VqnPKDNzKy6uInezMysCjnBm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqkBO8mZlZFXKCNzMzq0JO8GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7wZmZmVUgRUe4YikLSK8ALZXjpnYFXy/C61cSfYff5M+w+f4bd58+wOLryOe4dEcNybaiaBF8ukhZERGO54+jN/Bl2nz/D7vNn2H3+DIujWJ+jm+jNzMyqkBO8mZlZFXKC774flDuAKuDPsPv8GXafP8Pu82dYHEX5HN0Hb2ZmVoVcgzczM6tCTvAFknS8pKWSlkm6uNzx9DaS9pT0G0nPSlos6Qvljqm3klQjaZGkX5Y7lt5KUr2kWZL+mP5OHlHumHobSRek/5efkfQTSQPLHVOlk/QjSS9Leiaj7B2Sfi3pufTfoYWe3wm+AJJqgBuADwMHAadJOqi8UfU6m4CLIuJA4H3Aef4MC/YF4NlyB9HL/Qfwq4g4ADgUf55dIqkB+DzQGBHvAmqAU8sbVa9wK3B8VtnFwAMRsR/wQPq8IE7whTkMWBYRz0fEW8BMYHyZY+pVIuKliHg8ffwGyR/UhvJG1ftI2gP4Z+DmcsfSW0naATga+CFARLwVEc3ljapX6g/USeoPDAJWlzmeihcRvwVeyyoeD8xIH88AJhR6fif4wjQAKzOer8LJqWCSRgBjgEfLG0mvdB3w78DmcgfSi70TeAW4Je3quFnS4HIH1ZtERBNwLfAi8BLwekTcV96oeq1dI+IlSCpCwC6FnsgJvjDKUebbEQogaQhwF/DFiFhb7nh6E0knAC9HxMJyx9LL9QfeA9wUEWOAN+lGs2hflPYTjwdGAsOBwZLOLG9U5gRfmFXAnhnP98DNUV0mqZYkuf84Iu4udzy90JHAiZJWkHQTHSvp9vKG1CutAlZFRFsL0iyShG/5+xDw54h4JSI2AncD/1DmmHqrv0raHSD99+VCT+QEX5jHgP0kjZQ0gGQwyZwyx9SrSBJJn+ezEfHtcsfTG0XElIjYIyJGkPwOPhgRrjV1UUT8BVgpaVRadBywpIwh9UYvAu+TNCj9v30cHqhYqDnApPTxJOB/Cj1R/6KE08dExCZJ5wNzSUaL/igiFpc5rN7mSOCTwNOSnkjLvhIR95YxJuu7Pgf8OL1gfx74lzLH06tExKOSZgGPk9whswjPatcpST8BjgF2lrQKuBy4GrhT0tkkF04TCz6/Z7IzMzOrPm6iNzMzq0JO8GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb1YCkkLStzKef0nS1B6O4VZJp6SPb+7uYj6SRmSuelUskqZJ+lCO8mO6s0KepBWSdm5nm9J/p7Y9z1WW/vvjdOXIZ9LVv2oLjcmsJznBm5XG34GT2kswnUkX7CiaiPh0RJR88pZ0pcUuiYjLIuL+UsTTgXdL+i7wDkkTgK+1UwbwY+AA4GCgDvh0D8dqVhBPdGNWGptIJvq4APhq5gZJewM/AoaRLHLyLxHxoqRbSVaWGgM8LukNkrm9dwf2By4kWVr3w0AT8NGI2CjpMuCjJMnnd8BnImuCC0nzgC+RzBM+LS2uAwZExEhJ7wW+DQwBXgU+FREvpeU/AtYD/5frjUo6hmSCjpeAdwMHpfOQfx4YQLKI0Lnp7j8EGknWbvhRRHwnfd+/jIhZko4nWUDnVZJJU9peYyqwLiKuTZ8/A5wQESskzSaZOnog8B8RsdUEK+nCMXeSTCldA1wZET+V1AI8AtRGxL+l+25Tljn5kqQ/pOcxq3iuwZuVzg3AGZJ2zCq/HrgtIg4hqR1+N2Pb/sCHIuKi9Pk+JMvBjgduB34TEQcDLWk5wPURMTZdh7sOOKG9gCJiTkS8OyLeDTwJXJs2OX8POCUi2hJ6W+31FuDzEXFEJ+/1MOCrEXGQpAOBTwBHpq/TCpxBkvwbIuJd6Xu4JfMEkgYC/0VysfJ+YLdOXrPNv6ZxNwKfl7RT1vbjgdURcWj6Gf1K0rtJLjpuB+ZKuipXWVZ8tSSzL/4qz7jMysoJ3qxE0tXxbiOpyWY6ArgjffzfwFEZ234WEa0Zz/83XbzjaZLaZ1tyeRoYkT7+oKRHJT0NHAuM7iw2Sf8OtETEDcAo4F3Ar9Npgy8B9kgvTOojYn5GrO35Q0T8OX18HPBe4LH0fMeRLMn6PPBOSd9La+rZqwceQLJgyXNpC0S+C+d8XtKTwO9JavL7ZW1/GviQpGskvT8iXgeejIjPA3+LiNnApe2UZboR+G1EPJRnXGZl5SZ6s9K6jqSp+ZYO9slsTn8za9vfASJis6SNGU3vm4H+aa33RqAxIlamTdkDOwpI0nEk81sf3VYELM6upUuqJ/9lkDPjFjAjIqbkeO1DgXHAecDHgX/N2qW919vE1hWSgen5jiFZyeyIiFifdkVs9f4j4k9pV8NHgG9Iui8ipqXbpqb/Rsb+25RJupykS+Uz7cRnVnFcgzcroYh4jaT/9+yM4t+RrP4GSdN1zr7tPLUls1clDQFO6WjntP//RuDjEdGSFi8Fhkk6It2nVtLoiGgGXpfU1sJwRp4xPQCcImmX9HzvkLR3OuCwX0TcRVI7zl6S9Y/ASEn7pM9Py9i2om1/Se8hGZsAsCOwJk3uB5CMUch+z8OB9RFxO3BtjtftkKRPk1yUnBYRm7tyrFk5uQZvVnrfAs7PeP554EeSJpMOsiv0xBHRLOm/SJqhV5AsZdyRTwE7AT9P7wJbHREfSW+n+27aLN+fpOVhcRrbjyStJ1k9MZ+Ylki6BLhPUj9gI0mNvQW4JS0DmJJ13AZJ5wD3SHqV5MLnXenmu4Cz0ib/x4A/peW/Aj4r6SmSC5Xf5wjpYGC6pM1pLP+Wz/vI8J/AC8Aj6Wd2d1sLgFkl82pyZmZmVchN9GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQv8fHB2xvJAY948AAAAASUVORK5CYII=\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.1792145490062"
      ]
     },
     "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.893213374093957),\n",
       " ('p-value', 0.08658690502352209),\n",
       " ('f-value', 2.503715946256434),\n",
       " ('f p-value', 0.08794028782673029)]"
      ]
     },
     "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.1002422436378152), ('p-value', 0.3820295068692507)]"
      ]
     },
     "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.0796490077823977), ('p value', 0.28346392475408044)]"
      ]
     },
     "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}