{
 "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:                Mon, 24 Feb 2020   Prob (F-statistic):           1.90e-08\n",
      "Time:                        22:49:06   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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5zUdd3//8fTZZGFRRcTT4sCeQAlD9R6ujQz7fqhZUGapWlal369vNTUMlLKAx66tK9WVipWVsrXAyl6EZVXmKloZSaIJyAu0VBZLPDSFZElObx+f3w+i8Mwuzs7O7MzO/u83257Y+Zzmtd8ZtnXvM+KCMzMzKy6bFHuAMzMzKz4nODNzMyqkBO8mZlZFXKCNzMzq0JO8GZmZlXICd7MzKwKOcGbWYckTZZ0e/p4F0mrJNUU+TWWSPpYMa/ZzuvMl3R4O/sOl7S0SK/ziKTTi3GtSifpi5L+UO44bHNO8NYtPfWHuZql9/AfkgZlbDtd0iNlDCuniHglIuojYn25YylERIyJiEfKHYdZT3CCt15NUr9yx1Ak/YDzunsRJar6/3Wxaw/sPVX0/8lwgrcSknSMpKcltUj6k6R90u0XSZqedez3Jf0gfby1pJ9Kek1Ss6Sr2v6op9WBf5T0PUlvAJMl7SrpIUn/K+l1SXdIasi49gclzZP0tqR7JP1C0lWdxZnj/dws6bqsbb+U9NX08YVpvG9LWiTpyC7crmuBr2XGnfU6/yLpSUlvpf/+S8a+RyR9S9IfgdXA+9NtV6XvZ5WkX0l6X3pvVqbXGJF1/19N982V9OF24hghKST1k3Rweu22nzWSlqTHbZF+zi+mn8vdkrbJuM4XJL2c7vtmRzdG0q2Spki6X9I7wEclbSnpOkmvpLUfN0uqS4/fVtKv08/zDUmPtX3pyaxxklSXXvtNSQuA/bNeNyTtlhXHVenjIelrrEjP/7WkYe3Ev5uk2eln97qkX7Rz3ABJt6f3pCX9jLZP941Mr/G2pN9JukHvNZts1rSQ9T4PkPR4es3X0nP7Z73PsyW9ALyQbhudvs4b6e/yZzOOf5+kmenvyl+AXTv6/Kx8nOCtJCR9EPgZ8O/A+4AfATMlbQncBXxc0lbpsTXAZ4E709NvA9YBuwFjgf8PyGzPPBB4CdgO+BYg4GpgJ2BPYGdgcnrt/sB/AbcC26Sv/ek848x2J/A5SUrPHZLGNk3SKOAcYP+IGAyMA5Z04ZbNAR4Bvpa9I02MvwF+kMb4XeA3kt6XcdgXgDOAwcDL6bYT0u2NJH+EHwd+nt6HhcBlGec/CeyX7rsTuEfSgI4CjojH0+r6emAI8GeS+wtwLjAB+AjJ5/ImcGP6fvYCpqSx7ZS+p5zJMcPnST7rwcAfgG8De6Qx75a+x0vTYy8AlgJDge2BbwC55uS+jOS+7EryeZ3aSQyZtiC5l8OBXYBW4IZ2jr0SeIDkHg0DftjOcacCW5P8/r4PODO9LiSfyVxg2/R6XYl1PfCV9NyDgSOBs7KOmUDy/2ovJU1Fv0tfczvgROAmSWPSY28E1gA7Av+W/lgligj/+KfgH5Ik9rEc26cAV2ZtWwR8JH38B+CU9PG/Ai+mj7cH/gnUZZx3IvBw+viLwCudxDQBmJc+PgxoBpSx/w/AVfnEmbVdwCvAYenz/wM8lD7eDVgOfAyoLeQeAh8A3iJJTKcDj6T7vwD8Jeucx4Evpo8fAa7I2v8I8M2M598B/jvj+SeBpzuI6U1g3/TxZOD29PEIkmTZL8fn/Rtgi/T5QuDIjP07AmtJmiIuBaZl7BsEvJvr9yjdfyswNetzeAfYNWPbwcDf0sdXAL8Eduvo95XkS+JRGfvOAJZmPI/Ma6RxXNVOjPsBb2bd/9PTx1OBHwPDOvk9+DfgT8A+Wdt3IfnCOyhj250Zn8nhmXF39P8y3Xc+8F9Z7/OIjOefAx7LOudHJF+IatLPcXTGvv8E/tCV33n/9MyPS/BWKsOBC9JqwRZJLSQlk53S/XeSJG5ISmd3ZpxXC7yWcd6PSEoSbV7NfCFJ20mapqR6fCVwO0lphfT1miP9S5Tj/M7i3Ci9xrSsuO9I9y0m+cM5GViexrPZNToSEc8DvwYuytq1E++Vytu8TFJqzfWe2vwj43Frjuf1bU8kXSBpYVqN3EJSktyWPEj6d5Ik8/mI2JBuHg78V8Y9XUhSktw+fT8b442Id4D/7eRlMt/fUGAgMDfj+r9Nt0PS3LEYeEDSS5Ky72ebTeJg83vcLkkDJf0obWZYCTwKNCh3/4Cvk3wp+YuSXvztlXj/HzCLpEZomaT/K6k2jfPN9D4VEuseaRPC39NY/5PNP9vs/xMHZv2fOAnYgeQe96PA+2Y9ywneSuVV4FsR0ZDxMzAi2qpw7wEOT9stP817Cf5VkhL8thnnbRURYzKunV3denW6bZ+I2Ao4meQPKsBrQGNbtXpq5y7Eme0u4DOShpNUad67MaiIOyPiUJI/kEFSjdxVl5HUDGQm72XpNTPtQlIzsfHlC3gtAJS0t19I0kwyJCIaSGoS1OGJ7517JTA+It7K2PUqcHTWfR0QEc0kn8nOGdcYSFIl3ZHM9/c6yReUMRnX3jqSpgIi4u2IuCAi3k9SU/FV5e4PsUkcJPc002qSLxJtdsh4fAEwCjgw/Z07rO3tbBZ4xN8j4v9ExE4kTUE3ZbbtZxy3NiIuj4i9gH8BjgFOSeMcooxRFlmxvpMZZ/olY2jG/inAX4Hd01i/kSPO7C/As7M+u/qI+A9gBUltQkf3zSqEE7wVQ23aQajtpx/wE+BMSQcqMUjSJyQNBoiIFSTVmD8nqVpdmG5/jaS98juStlLSWWtXSR/p4PUHA6uAFkmNwMSMfY+TlBzPUdIxbDxwQMb+DuPMFhHzSP7I3QLMiogWAEmjJB2Rtt2vIUlAXR5KltYE/IKkDbvN/cAekj6fvofPAXuRlPaLYTDJH+0VQD9JlwJbdXaSpJ3TWE+JiP/J2n0z8K30ixCShqb3HmA6cIykQ9M+ElfQhb9FaS3BT4DvSdouvX6jpHHp42OUdGwTsJLkc8j1WdwNTFLSYW4Y8OWs/U8Dn5dUI+kokv4EbQaTfMYtaR+Jy2iHpOP1Xge8N0mS6WbxSPqopL3TBL2SpCp8fUS8TNJH43JJ/SUdSvLFpc3/AAPS39ta4GIgsw/J4PR6qySNBv6jvVhTvyb5ffuCpNr0Z39Je0YyPPI+ks6tA9P+FF3pD2A9yAneiuF+kj92bT+TI2IOSUn0BpI/aotJ2s8z3UnS9nxn1vZTgP7AgvTc6SRtuO25HPggSanzNyR/gACIiHeBY4HTgBaS0v2vSWoJyDPObHfliHtL4BqS0uXfSZoUvgEg6SRJ8zu5ZqYrSNql297D/5KU5i4gqcr+OnBMRLzehWt2ZBbw3ySJ4mWSLyi5qvyzHUlSqp2u93rSt73P7wMzSarJ3ybpgHdg+n7mA2eT3L/XSO57VyeYuZDks/pzWu38IEmJGmD39Pkqki94N0Xuse+Xk7zfv5F8qfx/WfvPI0mkbVXUMzL2XQ/UkXzefyZpImjP/sATklaR3JPzIuJvOY7bgeR3fSVJk8ZskuYmSJqDDgTeIPkyMbXtpLTm5CySL53NJCX6zPv5tfT8t0m+GOXsxZ9xvbdJOo+eQFJ79HeS2qi2Lw3nkDTv/J2kX8LPO7qelY82bZo0q36SngBujgj/YbJeSdJkkg6AJ5c7FqtcLsFb1ZP0EUk7pNXbpwL70HGJy8ys1/OsRdYXjCJpb60HXgQ+k7b1m5lVLVfRm5mZVSFX0ZuZmVUhJ3gzM7MqVFVt8Ntuu22MGDGi3GGYmZn1iLlz574eEUNz7auqBD9ixAjmzJlT7jDMzMx6hKR2pwp2Fb2ZmVkVcoI3MzOrQk7wVeKFF15gwIABnHyyJ7YyMzMn+Kpx9tlns//++5c7DDMzqxBO8FVg2rRpNDQ0cOSRuVbENDOzvsgJvpdbuXIll156Kd/5znfKHYqZmVUQJ/he7pJLLuG0005j5513LncoZmZWQapqHHxf8/TTT/Pggw8yb968codiZmYVxgm+F3vkkUdYsmQJu+yyCwCrVq1i/fr1LFiwgKeeeqrM0ZmZWTlV1WpyTU1N0Zdmslu9ejUrV67c+Py6665jyZIlTJkyhaFDc85caGZmVUTS3IhoyrXPJfhebODAgQwcOHDj8/r6egYMGODkbmZmTvDVZPLkyeUOwczMKoR70ZuZmVUhJ3gzM7MqVNIEL+koSYskLZZ0UY79oyU9Lumfkr6WY3+NpHmSfl3KOM3MzKpNyRK8pBrgRuBoYC/gREl7ZR32BnAucF07lzkPWFiqGM3MzKpVKUvwBwCLI+KliHgXmAaMzzwgIpZHxJPA2uyTJQ0DPgHcUsIYe70Z85o55JqHGHnRbzjkmoeYMa+53CGZmVkFKGWCbwRezXi+NN2Wr+uBrwMbOjpI0hmS5kias2LFiq5H2YvNmNfMpPueo7mllQCaW1qZdN9zTvJmZlbSBK8c2/KaVUfSMcDyiJjb2bER8eOIaIqIpr42/vvaWYtoXbt+k22ta9dz7axFZYrIzMwqRSkT/FIgcwWUYcCyPM89BPiUpCUkVftHSLq9uOH1fstaWru03czM+o5SJvgngd0ljZTUHzgBmJnPiRExKSKGRcSI9LyHIuLk0oXaO+3UUNel7WZm1neULMFHxDrgHGAWSU/4uyNivqQzJZ0JIGkHSUuBrwIXS1oqaatSxVRtJo4bRV1tzSbb6mprmDhuVJkiMjOzSuHFZnq5GfOauXbWIpa1tLJTQx0Tx41iwtiu9GU0M7PeyovNVLEJYxud0M3MbDOeqtbMzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7weaqvr9/kp6amhi9/+csb9//+979n9OjRDBw4kI9+9KO8/PLLZYzWzMz6Oif4PK1atWrjzz/+8Q/q6uo4/vjjAXj99dc59thjufLKK3njjTdoamric5/7XJkjNjOzvswJvgDTp09nu+2248Mf/jAA9913H2PGjOH4449nwIABTJ48mWeeeYa//vWvZY7UzMz6Kif4Atx2222ccsopSMmS9/Pnz2fffffduH/QoEHsuuuuzJ8/v1whmplZH+cE30WvvPIKs2fP5tRTT924bdWqVWy99dabHLf11lvz9ttv93R4ZmZmgBN8l02dOpVDDz2UkSNHbtxWX1/PypUrNzlu5cqVDB48uKfDMzMzA5zgu2zq1KmblN4BxowZwzPPPLPx+TvvvMOLL77ImDFjejo8MzMzwAm+S/70pz/R3Ny8sfd8m09/+tM8//zz3HvvvaxZs4YrrriCffbZh9GjR5cpUjMz6+uc4Lvgtttu49hjj92s6n3o0KHce++9fPOb32TIkCE88cQTTJs2rUxRmpmZgSKi3DEUTVNTU8yZM6fcYZiZmfUISXMjoinXPpfgzczMqpATvJmZWRVygjczM6tCTvBmZmZVqF+5A+jNZsxr5tpZi1jW0spODXVMHDeKCWMbyx2WmZmZE3yhZsxrZtJ9z9G6dj0AzS2tTLrvOQAneTMzKztX0Rfo2lmLNib3Nq1r13PtrEVlisjMzOw9TvAFWtbS2qXtZmZmPckJvkA7NdR1abuZmVlPcoIv0MRxo6irrdlkW11tDRPHjSpTRGZmZu9xJ7sCtXWkcy96MzOrRE7w3TBhbKMTupmZVSRX0ZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRUqaYKXdJSkRZIWS7oox/7Rkh6X9E9JX8vYvrOkhyUtlDRf0nmljNPMzKzalGyYnKQa4EbgX4GlwJOSZkbEgozD3gDOBSZknb4OuCAinpI0GJgr6XdZ55qZmVk7SlmCPwBYHBEvRcS7wDRgfOYBEbE8Ip4E1mZtfy0inkofvw0sBDzg3MzMLE+lTPCNwKsZz5dSQJKWNAIYCzxRlKjMzMz6gFImeOXYFl26gFQP3AucHxEr2znmDElzJM1ZsWJFAWGamZlVn1Im+KXAzhnPhwHL8j1ZUi1Jcr8jIu5r77iI+HFENEVE09ChQwsO1szMrJqUMsE/CewuaaSk/sAJwMx8TpQk4KfAwoj4bgljNDMzq0ol60UfEesknQPMAmqAn0XEfElnpvtvlrQDMAfYCtgg6XxgL2Af4AvAc5KeTi/5jYi4v1TxmpmZVZOSriaXJuT7s7bdnPH47yRV99n+QO42fDMzM8uDZ7IzMzOrQk7wZmZmVcgJ3szMrAo5wZuZmVUhJ3gzM7Mq5ARvZmZWhZzgzczMqpATvJmZWRVygjczM6tCTvBmZmZVyAnezMysCjnBm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqkBO8mZlZFXKCNzMzq0JO8Hk6/PDDGTBgAPX19dTX1zNq1KiN+1avXs1ZZ53Ftttuy9Zbb81hhx1WxkjNzMygX7kD6E1uuOEGTj/99M22n3HGGaxbt46FCxeyzTbb8PTTT5chOjMzs/c4wXfTokWLmDlzJkuXLmWrrbYC4EMf+lCZozIzs77OVfRdMGnSJLbddlsOOeQQHnnkEQCeeOIJhg8fzmWXXca2227L3nvvzb333lveQM3MrM9zgs/Tt7/9bV566SWam5s544wz+OQnP8mLL77I0qVLef7559l6661ZtmwZN9xwA6eeeioLFy4sd8hmZtaHOcHn6cADD2Tw4MFsueWWnHrqqRxyyCHcf//91NXVUVtby8UXX0z//v35yEc+wkc/+lEeeOCBcodsZmZ9mBN8gSQREeyzzz7lDsXMzGwzTvB5aGlpYdasWaxZs4Z169Zxxx138OijjzJu3DgOO+wwdtllF66++mrWrVvHH//4Rx555BHGjRtX7rDNzKwPcy/6PKxdu5aLL76Yv/71r9TU1DB69GhmzJixcSz8L3/5S04//XSuueYahg8fztSpUxk9enSZozYzs75MEVHuGIqmqakp5syZU+4wzMzMeoSkuRHRlGufq+jNzMyqkBO8mZlZFXKCNzMzq0LuZFegGfOauXbWIpa1tLJTQx0Tx41iwtjGcodlZmYGOMEXZMa8Zibd9xyta9cD0NzSyqT7ngNwkjczs4rgKvoCXDtr0cbk3qZ17XqunbWoTBGZmZltygm+AMtaWru03czMrKc5wRdgp4a6Lm03MzPraU7wBZg4bhR1tTWbbKurrWHiuFFlisjMzGxTJU3wko6StEjSYkkX5dg/WtLjkv4p6WtdObecJoxt5Opj96axoQ4BjQ11XH3s3u5gZ2ZmFaNkvegl1QA3Av8KLAWelDQzIhZkHPYGcC4woYBzy2rC2EYndDMzq1ilLMEfACyOiJci4l1gGjA+84CIWB4RTwJru3qumZmZta+UCb4ReDXj+dJ0W1HPlXSGpDmS5qxYsaKgQM3MzKpNKRO8cmzLd+m6vM+NiB9HRFNENA0dOjTv4MzMzKpZKRP8UmDnjOfDgGU9cK6ZmVmfV8oE/ySwu6SRkvoDJwAze+BcMzOzPq9kvegjYp2kc4BZQA3ws4iYL+nMdP/NknYA5gBbARsknQ/sFRErc51bqljNzMyqjSLybRavfE1NTTFnzpxyh2FmZtYjJM2NiKZc+zyTnZmZWRVygjczM6tCTvBmZmZVyAm+CF544QUGDBjAySefDMCSJUuQRH19/cafK6+8ssxRmplZX1KyXvR9ydlnn83++++/2faWlhb69fMtNjOznucSfDdNmzaNhoYGjjzyyHKHYmZmtpETfDesXLmSSy+9lO985zs59w8fPpxhw4bxpS99iddff72HozMzs77MCb4bLrnkEk477TR23nnnTbZvu+22PPnkk7z88svMnTuXt99+m5NOOqlMUZqZWV/kBuICPf300zz44IPMmzdvs3319fU0NSXzDmy//fbccMMN7LjjjqxcuZKtttqqp0M1M7M+yAm+QI888ghLlixhl112AWDVqlWsX7+eBQsW8NRTT21yrJQsjldNswaamVll81S1BVq9ejUrV67c+Py6665jyZIlTJkyhZdeeomGhgZ233133nzzTc466yyWL1/Oww8/3COxmZlZ3+Cpaktg4MCB7LDDDht/6uvrGTBgAEOHDuWll17iqKOOYvDgwXzgAx9gyy235K677ip3yGZm1oe4BG9mZtZLuQRvZmbWxzjBm5mZVSEneDMzsyrkYXJFNGNeM9fOWsSyllZ2aqhj4rhRTBjbWO6wzMysD3KCL5IZ85qZdN9ztK5dD0BzSyuT7nsOwEnezMx6nKvoi+TaWYs2Jvc2rWvXc+2sRWWKyMzM+jIn+CJZ1tLape1mZmal5ARfJDs11HVpu5mZWSk5wRfJxHGjqKut2WRbXW0NE8eNKlNEZmbWl7mTXZG0daRzL3ozM6sETvBFNGFsoxO6mZlVBFfRm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqUF4JXtL2kn4q6b/T53tJOq20oZmZmVmh8i3B3wrMAnZKn/8PcH4pAjIzM7PuyzfBbxsRdwMbACJiHbC+41PMzMysXPJN8O9Ieh8QAJIOAt4qWVRmZmbWLflOdPNVYCawq6Q/AkOBz5QsKjMzM+uWvBJ8RDwl6SPAKEDAoohYW9LIzMzMrGB5JXhJx2Zt2kPSW8BzEbG8+GGZmZlZd+TbBn8acAtwUvrzE5Jq+z9K+kJ7J0k6StIiSYslXZRjvyT9IN3/rKQPZuz7iqT5kp6XdJekAV16Z2ZmZn1Yvgl+A7BnRBwXEccBewH/BA4ELsx1gqQa4Ebg6PT4EyXtlXXY0cDu6c8ZwJT03EbgXKApIj4A1AAndOF9mZmZ9Wn5JvgREfGPjOfLgT0i4g2gvbb4A4DFEfFSRLwLTAPGZx0zHpgaiT8DDZJ2TPf1A+ok9QMGAsvyjNXMzKzPy7cX/WOSfg3ckz4/DnhU0iCgpZ1zGoFXM54vJSnxd3ZMY0TMkXQd8ArQCjwQEQ/kGauZmVmfl28J/myS2ez2A8YCU4GzI+KdiPhoO+cox7bI5xhJQ0hK9yNJZs8bJOnknC8inSFpjqQ5K1as6PydmJmZ9QF5Jfi0Cn16RHwlIs5PH2cn62xLgZ0zng9j82r29o75GPC3iFiRDse7D/iXdmL7cUQ0RUTT0KFD83k7Fenkk09mxx13ZKuttmKPPfbglltuKXdIZmbWi+W72MxBkp6UtErSu5LWS1rZyWlPArtLGimpP0knuZlZx8wETkl70x8EvBURr5FUzR8kaaAkAUcCC7v0znqZSZMmsWTJElauXMnMmTO5+OKLmTt3brnDMjOzXirfKvobgBOBF4A64HTghx2dkM5Xfw7JIjULgbsjYr6kMyWdmR52P/ASsJhk6N1Z6blPANOBp4Dn0jh/nP/b6n3GjBnDlltuCYAkJPHiiy+WOSozM+ut1HlNO0iaExFNkp6NiH3SbX+KiJzV5uXS1NQUc+bMKXcYBTvrrLO49dZbaW1tZezYsTz66KPU19eXOywzM6tQkuZGRFOuffmW4Fen1exPS/q/kr4CDCpahAbATTfdxNtvv81jjz3Gscceu7FEb2Zm1lX5JvgvpMeeA7xD0jHuuFIF1ZfV1NRw6KGHsnTpUqZMmVLucMzMLIdp06ax5557MmjQIHbddVcee+wxFixYQFNTE0OGDGHIkCF87GMfY8GCBWWLsdNx8OmMdN+KiJOBNcDlJY/KWLdundvgzcwq0O9+9zsuvPBCfvGLX3DAAQfw2muvATBo0CCmT5/O8OHD2bBhAzfeeCMnnHACzz77bFni7LQEHxHrgaFpFb2VwPLly5k2bRqrVq1i/fr1zJo1i7vuuosjjjii3KGZmVmWyy67jEsvvZSDDjqILbbYgsbGRhobG2loaGDEiBFIIiKoqalh8eLFZYsz35nslpAsLDOTpIoegIj4bimC6mskMWXKFM4880w2bNjA8OHDuf766xk/PntmXzMzK6f169czZ84cPvWpT7HbbruxZs0aJkyYwLXXXktdXR0ADQ0NrFq1ig0bNnDFFVeULdZ8E/yy9GcLYHDpwumbhg4dyuzZs8sdhpmZdeIf//gHa9euZfr06Tz22GPU1tYyfvx4rrrqKr71rW8B0NLSwjvvvMNtt93G8OHDyxZrXsPkNh4sDYqIdzo/sjx6+zA5MzOrbG+++SbbbLMNt956K6eeeioA9957L1dddRXz5s3b5NgNGzYwdOhQFi5cyHbbbVeSeLo9TE7SwZIWkM4mJ2lfSTcVMUYzM7OKN2TIEIYNG0YyyWrHNmzYwOrVq2lubu6ByDaX7zC564FxwP8CRMQzwGGlCsrMzKxSfelLX+KHP/why5cv58033+T666/nmGOO4Xe/+x3z5s1j/fr1rFy5kq9+9asMGTKEPffcsyxx5pvgiYhXszatL3IsfdqMec0ccs1DjLzoNxxyzUPMmFeeb3xmZtaxSy65hP3335899tiDPffck7Fjx/LNb36TlpYWTjzxRLbeemt23XVXFi9ezG9/+1sGDBhQljjznap2OvBdkjnpDwLOBZoi4oTShtc1vbUNfsa8Zibd9xyta9/7zlRXW8PVx+7NhLGNZYzMzMwqWTGmqj2TZE34RpIlXvdLn1sRXDtr0SbJHaB17XqunbWoTBGZmVlvl+8wOUXESSWNpA9b1tLape1mZmadybcE/ydJD0g6TVJDSSPqg3ZqqOvSdjMzs87kleAjYnfgYmAM8JSkX0s6uaSR9SETx42irrZmk211tTVMHDeqTBGZmVlvl28VPRHxF+Avkv6TpMPdbcDtpQqsL2nrSHftrEUsa2llp4Y6Jo4b5Q52Zma9wIx5zRX59zuvBC9pK+DTwAnArsB/AQeUMK4+Z8LYxor4hTAzs/xlj4Jqbmll0n3PAZT9b3q+bfDPkPScvyIi9oiICyNibgnjMjMzq3iVPAoq3yr690dESBpU0mjMzMx6kUoeBZVvCf4gz0VvZma2qUoeBeW56M3MzApUyaOgutKL/tWs1XM8F72ZmfVplTwKKt8E/6qkfwFCUn+SuegXli4s60ilDskwM+uLKnUUVL4J/kzg+7w3F/0DwFmlCsraV8lDMszMrHLkO5Pd6xFxUhaVHtMAAB3ASURBVERsHxHbRcTJwCkljs1yqOQhGWZmVjnyXg8+h68WLQrLW3tDL5orYEiGmZlVju4keHV+iBVbe0MvRFJ9b2ZmBt1L8FG0KCxvE8eNyvnNKsDV9GZmtlGHCV7S25JW5vh5G9iph2K0DBPGNrb7zaoSZk4yM7PK0GEv+ogY3FOBWP4aG+pytrlXwsxJZmZWGbpTRW9lUskzJ5mZWWXIeyY7qxyVPHOSmZlVBif4XqpSZ04yM7PK4Cp6MzOzKuQSvPUoz6NvZtYznOCtx3gefTOznlPSKnpJR0laJGmxpIty7JekH6T7n5X0wYx9DZKmS/qrpIWSDi5lrFZ6nkffzKznlCzBS6oBbgSOBvYCTpS0V9ZhRwO7pz9nAFMy9n0f+G1EjAb2xcvT9nrtTcTjCXrMzIqvlFX0BwCLI+IlAEnTgPHAgoxjxgNTIyKAP6el9h2Bd4DDgC8CRMS7wLsljNV6wE4VOEGP+wSYWbUqZRV9I/BqxvOl6bZ8jnk/sAL4uaR5km6RNCjXi0g6Q9IcSXNWrFhRvOit6Cptgp62PgHNLa0E7/UJ8KI9ZlYNSpng21sTJZ9j+gEfBKZExFiSEv1mbfgAEfHjiGiKiKahQ4d2J96imDGvmUOueYiRF/2GQ655yMkiw4SxjVx97N40NtQhkil3rz5277KVmN0nwMyqWSmr6JcCO2c8HwYsy/OYAJZGxBPp9um0k+AriXuJd66SJuhxnwAzq2alLME/CewuaaSk/sAJwMysY2YCp6S96Q8C3oqI1yLi78Crktrqbo9k07b7iuQSYe/SXtu/F+0xs2pQsgQfEeuAc4BZJD3g746I+ZLOlHRmetj9wEvAYuAnwFkZl/gycIekZ4H9gP8sVazF4hJh71JpfQLMzIqppBPdRMT9JEk8c9vNGY8DOLudc58GmkoZX7FVYi9xa58X7TGzauaZ7Ipo4rhRm7TBg0uEla6S+gSYmRWTE3wRuURoZmaVwgm+yFwiNDOzSuDlYs3MzKqQE7yZmVkVcoI3MzOrQk7wZmZmVcid7MrAK5iZmVmpOcH3MM9Xb2ZmPcFV9D3M89WbmVlPcILvYZ6v3szMeoITfA/zCmZmZtYTnOB70Ix5zax+d91m2z1fvZmZFZs72fWQ7M51bRrqapn8qTHuYGdmZkXlBF9EHQ1/y9W5DmDQlv02Se4eQmdmZsXgBF8knQ1/y6dznYfQmZlZsbgNvkg6G/6WT+c6D6EzM7NicYIvks5K6BPHjaKutmaTfdmd6zyEzszMisUJvkg6K6FPGNvI1cfuTWNDHQIaG+q4+ti9N6l69xA6MzMrFrfBF8nEcaM26yWfXUKfMLaxw7b0fK5hZmaWD5fgi2TC2EaO+1AjNRIANRLHfajjhJ7rGp2V8s3MzPLhEnyRzJjXzL1zm1kfAcD6CO6d20zT8G26nOR7U0L3sD4zs8rkBN9NbQmuOUdHuLYe8JWe8ApN0h7WZ2ZWuVxF3w1tCS5Xcm9T6T3gM99D8F6SnjGvudNzPazPzKxyuQTfDZNnzs85O12mfHvAl6uqu6Mk3dnre1ifmVnlcgm+QDPmNdPSurbDY/LtAd+dUnR3dSdJe1ifmVnlcoIvUGfV0F3pAV/Oqu7uJOl8Ju8xM7PycBV9gToq4V7/uf26VL1ezqru7oy9z1xIx73ozcwqixN8gXZqqMvZuW7IwNouJ7j2rtUTVd3dTdK9bVifmVlf4QRfoPZKvpd9ckzRrtVTVd1O0mZm1ccJvkDFrJ7uyapuT0xjZtY3KNKZ16pBU1NTzJkzp9xhVKzsiWkgqSnwdLhmZr2TpLkR0ZRrn0vwRVbJJeTujHk3M7PexQm+iHpq6tZCv0R4Yhozs77DCb4bshPt6nfXlbyE3J0vEeXsrW9mZj2rpBPdSDpK0iJJiyVdlGO/JP0g3f+spA9m7a+RNE/Sr0sZZyFyzT735urcM9sVs4TcnUlxPDGNmVnfUbIEL6kGuBE4GtgLOFHSXlmHHQ3snv6cAUzJ2n8esLBUMXZHrkTbnmKWkLtTze715s3M+o5SVtEfACyOiJcAJE0DxgMLMo4ZD0yNpCv/nyU1SNoxIl6TNAz4BPAt4KsljLMg+ZbKi1lCnjGvmS2kjWvOZ8r3S4THvJuZ9Q2lTPCNwKsZz5cCB+ZxTCPwGnA98HVgcAljLFh77dltlB5TrF70bU0CuZK7q9nzU8kjHMzMiq2UCV45tmVnp5zHSDoGWB4RcyUd3uGLSGeQVO+zyy67FBJnQSaOG8VXfvH0Zm8IkqrvP150RFFfr70mgRrJ1ex56KkRDmZmlaKUneyWAjtnPB8GLMvzmEOAT0laAkwDjpB0e64XiYgfR0RTRDQNHTq0WLF3asLYRk46aJfNvqGUqjTdXpPAhggnqDyUc8U+M7NyKGWCfxLYXdJISf2BE4CZWcfMBE5Je9MfBLwVEa9FxKSIGBYRI9LzHoqIk0sYa0GumrA33/vcfp12Wpsxr5lDrnmIkRf9hkOueaigdd699nr3eA4AM+trSlZFHxHrJJ0DzAJqgJ9FxHxJZ6b7bwbuBz4OLAZWA18qVTzFlt2e+712logtVtVwuRek6e08B4CZ9TWei74AXZnT/ZBrHsqZWAppp3cnscJ5Hn4zq0aei77IujKnezGrhj3ErXA9uWKfmVklcIIvQFeStquGK4e/IJlZX1LSqWqrVVc6vPWl6WGL0ZnQzMyKwwm+AF1J2n1lethcc/NPuu85J3kzszJxFX0Butqe2xeqhnvjWvPutGhm1cwJvkCFJO1qTii9bZy5Z7Yzs2rnKvoeMmNeMxOnP7NJFfbE6c9UTRV2b5uIxzPbmVm1c4LvIZf/aj5r128658Da9cHlv5pfpoiKq7d1JuxtNQ5mZl3lKvoe8ubqtZ1u781V+L1tnLmHL5pZtXOCrxDV0CbcmzoTeupfM6t2rqLvIQ11tR1ud5twz+orwxfNrO9yCb6HTP7UGCbe8wxrN7zXDl+7hZj8qTGA24TLoTfVOJiZdZVL8AXq6qxtE8Y2cu3x+25SYrz2+H03Jpj22n4bBtZ6djizbqqvr9/kp6amhi9/+csb999yyy3stttu1NfXc9RRR7Fs2bIyRmtWHF5NrgClWJks1zVrawTBJqV+r4Bm1j3vvPMO22+/Pffffz+HHXYYs2fP5vjjj+fhhx9m991357zzzmPBggXMnj273KGadaqj1eRcgi9AKdrLc7UJD+rfb5PkXozXMevrpk+fznbbbceHP/xhAH71q19x/PHHM2bMGPr3788ll1zCo48+yosvvljmSM26x23wBWivXTzXsKuuyG4THnnRb/J6/d48vM6sp912222ccsopSAIgIsisyWx7/Pzzz7PrrruWJUazYnAJvgDttZcLitpGns/scOVY5MWrxllv9corrzB79mxOPfXUjds+/vGPc/fdd/Pss8/S2trKFVdcgSRWr15dxkjNus8JvgDtjZUOKGr1eT6zw/X08Dp/obDebOrUqRx66KGMHDly47YjjzySyy+/nOOOO47hw4czYsQIBg8ezLBhw8oYqVn3OcEXWWfD2rqSrPIZq92Tw+tmzGvmgrufqfovFFa9pk6duknpvc3ZZ5/NCy+8wPLlyznuuONYt24dH/jAB8oQoVnxuA2+AJNntj9//BYSIy/6zca2cHhv+taGgbWsWrNuY8e5fGar62ysdk9NudqWaNe3M+qiVOP1e+MytFaZ/vSnP9Hc3Mzxxx+/yfY1a9awePFixowZw6uvvsoZZ5zBeeedx5AhQ8oUqVlxuARfgJbW3PPKA6yPeG+1uHue2WQFuTdXry16r/ieWuQlV6LNVKo53NvruOgJgKyrbrvtNo499lgGDx68yfY1a9bw+c9/nvr6eg444AAOPvhgrrzyyjJFaVY8LsGXUHYyb09zSyuHXPNQQb3gC13kpas97ztKqKWaw33GvGZE0rchmxeFsa760Y9+lHN7Q0MDzz77bA9HY1Z6TvAFGDKwtt3V4QrVVlItZJGZrk65WsjCNu01BdRIJZt459pZi3Imd9F+R0czM0u4ir4Al31yTEmvX+rJbArped9eU8B3PrtvydrC26s1CHrPCnvFdPLJJ7Pjjjuy1VZbsccee3DLLbdsdszll1+OJB588MEyRGhmlcQJvgATxja2uzpcptotlEw3W4DuTprTkUJ63pdj9bX2quEb+2j1/KRJk1iyZAkrV65k5syZXHzxxcydO3fj/hdffJHp06ez4447ljFKM6sUTvAFmvypMXSUutsWk7n2M5suMJOvGhX2xSAf+Uygk8uEsY388aIj+Ns1n+CPFx1R8lJ0T3Ug7C3GjBnDlltuCYAkJG0yneo555zDt7/9bfr371+uEHsVz69g1c4JvkATxjbmbB+GpI24LQFmJ8V8k3x7w9GKobckTq/ZvrmzzjqLgQMHMnr0aHbccUc+/vGPA3DPPffQv3//jc+tY55fwfoCd7LrhsYCxqBPHDdqs1Xj2rt2d7XXU77Qnvfl4DXbN3XTTTfxwx/+kMcff5xHHnmELbfcklWrVvGNb3yDBx54oNzh9RqeX8H6Aif4bsiVrDsrCWcn1+zJb/K5Rj466ynvxNl71dTUcOihh3L77bczZcoUXn75Zb7whS9sMv2qdawnZ4A0Kxcn+G4otCScnVy7OiY9n+NdQql+69at48UXX2T27NksXbqUm266CYAVK1bw2c9+lgsvvJALL7ywzFFWpp6aAdKsnJzgu6kYJeGuXCPfMey9rYTiJW87tnz5ch566CGOOeYY6urqePDBB7nrrru48847ufTSS1m79r15Gfbff3+++93vcvTRR5cx4spWSO2bWW/jBN/L5FsyL0UJpVRJuJCJd/oaSUyZMoUzzzyTDRs2MHz4cK6//nrGjx+/2bE1NTUMGTKE+vr6MkTaO/SmfihmhXKC72XyLZkXu4RSyiTs5oTODR06lNmzZ+d17JIlS0obTJVwPxSrdh4mVyKlGmOb7xj2Yg8xK+W6872tOcHMrDdwCb4EcpV2z//F05z/i6dpqKtl8qfGFJxou1IyL2YJpZRJ2B2ezMyKzyX4EuhoadWW1rVMvOeZLpfo22oEvvKLp9my3xYMGVjbo5O/FDr7XT56y8Q7Zma9SUkTvKSjJC2StFjSRTn2S9IP0v3PSvpgun1nSQ9LWihpvqTzShlnsXVWql27IbpUtZ0961ZL61rWrN3A9z63X49MGQulTcKesc7MrPhKVkUvqQa4EfhXYCnwpKSZEbEg47Cjgd3TnwOBKem/64ALIuIpSYOBuZJ+l3VuxWqvyjlTe18CcvVUr4ROaKXudewOT13noYVm1pFStsEfACyOiJcAJE0DxgOZSXo8MDUiAvizpAZJO0bEa8BrABHxtqSFQGPWuRUrn+loc1Vtt9dTvb3r9HQnNCfhyuGhhWbWmVIm+Ebg1YznS0lK550d00ia3AEkjQDGAk/kehFJZwBnAOyyyy7dDLnrOipFTZ45n5bWtZudU7uFclZtt1dSr5FyLj7jTmh9VyXU6phZZStlgs+13ml2lurwGEn1wL3A+RGxMteLRMSPgR8DNDU1lW4Jthzyme99xrxmLv/VfN5cnST6jnrRt1ciXx9BXW2NZ92yjTy00Mw6U8oEvxTYOeP5MGBZvsdIqiVJ7ndExH0ljLNg+ZSiulKt3V7bfUNd7cZrAwwZWMtln8z9JcHtsn1D2+/Kyrm/4p3nf8+7K5YwaM+PsO9J3wDg3Xff5fOf/zxz5szh5Zdf5uGHH+bwww8vb9Bm1qNK2Yv+SWB3SSMl9QdOAGZmHTMTOCXtTX8Q8FZEvCZJwE+BhRHx3RLG2C3FLkXl6qleu4V45911m1T1r1m7Ief5XuO672j7XelX/z62Pvhz1O/9r9RkNf20rTi3ww47lDFSMyuXkpXgI2KdpHOAWUAN8LOImC/pzHT/zcD9wMeBxcBq4Evp6YcAXwCek/R0uu0bEXF/qeLtqhnzmpMGhhyNAoW2jefqqb763XUbq/fbtK5dz+SZ80u2gpxrASrfxt+VQf1Z1tLKli1L2L3+3Y3b+/fvz/nnnw8kc9ObWd9T0pns0oR8f9a2mzMeB3B2jvP+QO72+YowY14zF9zzDDn6vbXbgS5f2VX6Iy/6Tc7jWlrXMmNec9FXkKum3tnV/kUl83fl4osfZ+nSpWWOyMwqiWeyK8Dlv5rP+g25+/PVD+hX1CTSUW1A9mQ5xZhtrpRzzvckN1eYWV/nBF+A7CrzTC0d7CtER7UBuVaQ6+5sc9XSO7tavqiYmRXKCb7Iij02fcLYRoYMrM3rtYox5Wsp55zvSdXyRcXMrFBO8AVoG7aWS3NLa1GXhwW47JNj8i6ZTxjbyB8vOoK/XfOJguapr5aFX6rli0pn1q1bx5o1a1i/fj3r169nzZo1rFu3DoB//vOfrFmzBkiGza1Zs4bI1XHEzKqSE3wBJn9qDLVbtN8HMFd7b3fWh++pxVjaOqW1zZ5HCV+r1Lr6RaU7n085XXXVVdTV1XHNNddw++23U1dXx1VXXQXAqFGjqKuro7m5mXHjxlFXV8fLL79c5ojNrKeomr7RNzU1xZw5c3rktS6e8Rx3/PmVXKPkNmpsqOOPFx2xWc90SJJNJSXOSoyxu73g8z2/Et+7mVk+JM2NiKZc+0o6TK6aPfzXFR0md3ivvbc3zBteaTEWY7hevrMIVtp7NzMrBif4AuXTWautvbeSO3y1lXLbW9423xiLPea8J5NuJX8+ZmaFcht8gTrrrJXZ3lupHb4yx4q3J58YSzHmvCeTbqV+PmZm3eEEX6BcnbjaZHdMK7Rneqk7fuUqJXc1xvau090x5z2ZdKtl5AD03s6CZlZ8TvAFmjC2keM+1LjZfLptiSGzGrmQXvA9MRNbR6XhIQNr8+5kVorSdk8m3Z4apVBqnr3PzDK5Db4bcnW0a6+duCvLxkLPtEG3tzwttL9iXVeu053Sdq6Fd0o5l3xXP59K5M6CZpbJCb4bStlO3BNt0BPHjdpseFibriSGXNcpRmm7GpJuT3JnQTPL5Cr6bihlO3FPtEG3VU23J9/EUC1V3L2dOwuaWSYn+G4oZTtxT7VBTxjbSGMREkN3p8i17qumzoJm1n1O8N1QypJrT5aKnRiqg2tSzCyTp6o1oPgT1ZiZWel5qlrrlDu0mZlVF1fRm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqkBO8mZlZFXKCNzMzq0JO8GZmZlXICd7MzKwKOcGbmZlVISd4MzOzKuQEb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7wZmZmVcgJ3szMrAopIsodQ9FIWgG83MMvuy3weg+/ZrXxPew+38Pu8z0sDt/H7uvKPRweEUNz7aiqBF8OkuZERFO54+jNfA+7z/ew+3wPi8P3sfuKdQ9dRW9mZlaFnODNzMyqkBN89/243AFUAd/D7vM97D7fw+Lwfey+otxDt8GbmZlVIZfgzczMqpATfIEkHSVpkaTFki4qdzy9jaSdJT0saaGk+ZLOK3dMvZWkGknzJP263LH0VpIaJE2X9Nf0d/LgcsfU20j6Svp/+XlJd0kaUO6YegNJP5O0XNLzGdu2kfQ7SS+k/w4p5NpO8AWQVAPcCBwN7AWcKGmv8kbV66wDLoiIPYGDgLN9Dwt2HrCw3EH0ct8HfhsRo4F98f3sEkmNwLlAU0R8AKgBTihvVL3GrcBRWdsuAn4fEbsDv0+fd5kTfGEOABZHxEsR8S4wDRhf5ph6lYh4LSKeSh+/TfIHtbG8UfU+koYBnwBuKXcsvZWkrYDDgJ8CRMS7EdFS3qh6pX5AnaR+wEBgWZnj6RUi4lHgjazN44Hb0se3ARMKubYTfGEagVczni/FyalgkkYAY4EnyhtJr3Q98HVgQ7kD6cXeD6wAfp42ddwiaVC5g+pNIqIZuA54BXgNeCsiHihvVL3a9hHxGiSFIWC7Qi7iBF8Y5djm4QgFkFQP3AucHxEryx1PbyLpGGB5RMwtdyy9XD/gg8CUiBgLvEOBVaJ9VdpGPB4YCewEDJJ0cnmjMif4wiwFds54PgxXR3WZpFqS5H5HRNxX7nh6oUOAT0laQtJMdISk28sbUq+0FFgaEW01SNNJEr7l72PA3yJiRUSsBe4D/qXMMfVm/5C0I0D67/JCLuIEX5gngd0ljZTUn6Qzycwyx9SrSBJJm+fCiPhuuePpjSJiUkQMi4gRJL+DD0WES01dFBF/B16VNCrddCSwoIwh9UavAAdJGpj+3z4Sd1TsjpnAqenjU4FfFnKRfkULpw+JiHWSzgFmkfQW/VlEzC9zWL3NIcAXgOckPZ1u+0ZE3F/GmKzv+jJwR/qF/SXgS2WOp1eJiCckTQeeIhkhMw/PaJcXSXcBhwPbSloKXAZcA9wt6TSSL0/HF3Rtz2RnZmZWfVxFb2ZmVoWc4M3MzKqQE7yZmVkVcoI3MzOrQk7wZmZmVcgJ3qwEJIWk72Q8/5qkyT0cw62SPpM+vqW7i/lIGpG54lWxSLpC0sdybD+8OyvkSVoiadt29in9d3Lb81zb0n/vSFeOfD5d+au20JjMepITvFlp/BM4tr0E05l0wY6iiYjTI6Lkk7ekKy12SURcGhEPliKeDuwn6QfANpImAN9qZxvAHcBoYG+gDji9h2M1K4gnujErjXUkE318Bfhm5g5Jw4GfAUNJFjn5UkS8IulWklWlxgJPSXqbZG7vHYE9gK+SLK17NNAMfDIi1kq6FPgkSfL5E/DvkTXBhaRHgK+RzBN+Rbq5DugfESMlfQj4LlAPvA58MSJeS7f/DFgN/CHXG5V0OMnkHK8B+wF7pfOQnwv0J1lE6Kz08J8CTSRrN/wsIr6Xvu9fR8R0SUeRLKDzOsmkKW2vMRlYFRHXpc+fB46JiCWSZpBMHT0A+H5EbDLBSrpwzN0kU0rXAFdGxC8ktQKPA7UR8R/psZtty5x8SdJf0uuYVTyX4M1K50bgJElbZ22/AZgaEfuQlA5/kLFvD+BjEXFB+nxXkuVgxwO3Aw9HxN5Aa7od4IaI2D9dh7sOOKa9gCJiZkTsFxH7Ac8A16VVzj8EPhMRbQm9rfT6c+DciDi4k/d6APDNiNhL0p7A54BD0tdZD5xEkvwbI+ID6Xv4eeYFJA0AfkLyZeXDwA6dvGabf0vjbgLOlfS+rP1HAcsiYt/0Hv1W0n4kXzpuB2ZJuirXtqz4aklmX/xtnnGZlZUTvFmJpKvjTSUpyWY6GLgzffz/gEMz9t0TEesznv93unjHcySlz7bk8hwwIn38UUlPSHoOOAIY01lskr4OtEbEjcAo4APA79Jpgy8GhqVfTBoiYnZGrO35S0T8LX18JPAh4Mn0ekeSLMn6EvB+ST9MS+rZqweOJlmw5IW0BiLfhXPOlfQM8GeSkvzuWfufAz4m6duSPhwRbwHPRMS5wP9GxAzgkna2ZboJeDQiHsszLrOychW9WWldT1LV/PMOjsmsTn8na98/ASJig6S1GVXvG4B+aan3JqApIl5Nq7IHdBSQpCNJ5rY+rG0TMD+7lC6pgfyXQc6MW8BtETEpx2vvC4wDzgY+C/xb1iHtvd46Ni2QDEivdzjJSmYHR8TqtClik/cfEf+TNjV8HLha0gMRcUW6b3L6b2Qcv9k2SZeRNKn8ezvxmVUcl+DNSigi3iBp/z0tY/OfSFZ/g6TqOmfbdp7aktnrkuqBz3R0cNr+fxPw2YhoTTcvAoZKOjg9plbSmIhoAd6S1FbDcFKeMf0e+Iyk7dLrbSNpeNrhcIuIuJekdJy9JOtfgZGSdk2fn5ixb0nb8ZI+SNI3AWBr4M00uY8m6aOQ/Z53AlZHxO3AdTlet0OSTif5UnJiRGzoyrlm5eQSvFnpfQc4J+P5ucDPJE0k7WRX6IUjokXST0iqoZeQLGXckS8C7wP+Kx0FtiwiPp4Op/tBWi3fj6TmYX4a288krSZZPTGfmBZIuhh4QNIWwFqSEnsr8PN0G8CkrPPWSDoD+I2k10m++Hwg3X0vcEpa5f8k8D/p9t8CZ0p6luSLyp9zhLQ3cK2kDWks/5HP+8hwM/Ay8Hh6z+5rqwEwq2ReTc7MzKwKuYrezMysCjnBm5mZVSEneDMzsyrkBG9mZlaFnODNzMyqkBO8mZlZFXKCNzMzq0JO8GZmZlXo/wej7Ny0ydCmHQAAAABJRU5ErkJggg==\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.2rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
