vincentarelbundock/example_ols.ipynb

## example_ols.ipynb
{
 "metadata": {
  "name": "example_ols"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Ordinary least squares"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import statsmodels.api as sm\n",
      "import matplotlib\n",
      "import matplotlib.pyplot as plt\n",
      "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
      "np.random.seed(9876789)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## OLS estimation\n",
      "\n",
      "Artificial data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nsample = 100\n",
      "x = np.linspace(0, 10, 100)\n",
      "X = np.column_stack((x, x**2))\n",
      "beta = np.array([1, 0.1, 10])\n",
      "e = np.random.normal(size=nsample)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Our model needs an intercept so we add a column of 1s:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = sm.add_constant(X)\n",
      "y = np.dot(X, beta) + e"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/local/lib/python2.7/dist-packages/statsmodels-0.5.0-py2.7-linux-x86_64.egg/statsmodels/tools/tools.py:306: FutureWarning: The default of `prepend` will be changed to True in 0.5.0, use explicit prepend\n",
        "  FutureWarning)\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Inspect data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = sm.add_constant(X)\n",
      "y = np.dot(X, beta) + e"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Fit and summary:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "model = sm.OLS(y, X)\n",
      "results = model.fit()\n",
      "print results.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "                            OLS Regression Results                            \n",
        "==============================================================================\n",
        "Dep. Variable:                      y   R-squared:                       0.968\n",
        "Model:                            OLS   Adj. R-squared:                  0.968\n",
        "Method:                 Least Squares   F-statistic:                     1479.\n",
        "Date:                Sun, 26 Aug 2012   Prob (F-statistic):           2.19e-73\n",
        "Time:                        20:52:52   Log-Likelihood:                -146.51\n",
        "No. Observations:                 100   AIC:                             299.0\n",
        "Df Residuals:                      97   BIC:                             306.8\n",
        "Df Model:                           2                                         \n",
        "==============================================================================\n",
        "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
        "------------------------------------------------------------------------------\n",
        "x1             0.8598      0.145      5.948      0.000         0.573     1.147\n",
        "x2             0.1103      0.014      7.883      0.000         0.082     0.138\n",
        "const         10.3423      0.313     33.072      0.000         9.722    10.963\n",
        "==============================================================================\n",
        "Omnibus:                        2.042   Durbin-Watson:                   2.274\n",
        "Prob(Omnibus):                  0.360   Jarque-Bera (JB):                1.875\n",
        "Skew:                           0.234   Prob(JB):                        0.392\n",
        "Kurtosis:                       2.519   Cond. No.                         144.\n",
        "==============================================================================\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Quantities of interest can be extracted directly from the fitted model. Type ``dir(results)`` for a full list. Here are some examples:  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Parameters: ', results.params\n",
      "print 'R2: ', results.rsquared"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Parameters:  [  0.85975052   0.11025357  10.34233516]\n",
        "R2:  0.968242518536\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## OLS non-linear curve but linear in parameters\n",
      "\n",
      "We simulate artificial data with a non-linear relationship between x and y:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nsample = 50\n",
      "sig = 0.5\n",
      "x = np.linspace(0, 20, nsample)\n",
      "X = np.c_[x, np.sin(x), (x-5)**2, np.ones(nsample)]\n",
      "beta = [0.5, 0.5, -0.02, 5.]\n",
      "y_true = np.dot(X, beta)\n",
      "y = y_true + sig * np.random.normal(size=nsample)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Fit and summary:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "res = sm.OLS(y, X).fit()\n",
      "print res.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "                            OLS Regression Results                            \n",
        "==============================================================================\n",
        "Dep. Variable:                      y   R-squared:                       0.933\n",
        "Model:                            OLS   Adj. R-squared:                  0.928\n",
        "Method:                 Least Squares   F-statistic:                     211.8\n",
        "Date:                Sun, 26 Aug 2012   Prob (F-statistic):           6.30e-27\n",
        "Time:                        20:52:52   Log-Likelihood:                -34.438\n",
        "No. Observations:                  50   AIC:                             76.88\n",
        "Df Residuals:                      46   BIC:                             84.52\n",
        "Df Model:                           3                                         \n",
        "==============================================================================\n",
        "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
        "------------------------------------------------------------------------------\n",
        "x1             0.4687      0.026     17.751      0.000         0.416     0.522\n",
        "x2             0.4836      0.104      4.659      0.000         0.275     0.693\n",
        "x3            -0.0174      0.002     -7.507      0.000        -0.022    -0.013\n",
        "const          5.2058      0.171     30.405      0.000         4.861     5.550\n",
        "==============================================================================\n",
        "Omnibus:                        0.655   Durbin-Watson:                   2.896\n",
        "Prob(Omnibus):                  0.721   Jarque-Bera (JB):                0.360\n",
        "Skew:                           0.207   Prob(JB):                        0.835\n",
        "Kurtosis:                       3.026   Cond. No.                         221.\n",
        "==============================================================================\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Extract other quantities of interest:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Parameters: ', res.params\n",
      "print 'Standard errors: ', res.bse\n",
      "print 'Predicted values: ', res.predict()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Parameters:  [ 0.46872448  0.48360119 -0.01740479  5.20584496]\n",
        "Standard errors:  [ 0.02640602  0.10380518  0.00231847  0.17121765]\n",
        "Predicted values:  [  4.77072516   5.22213464   5.63620761   5.98658823   6.25643234\n",
        "   6.44117491   6.54928009   6.60085051   6.62432454   6.6518039\n",
        "   6.71377946   6.83412169   7.02615877   7.29048685   7.61487206\n",
        "   7.97626054   8.34456611   8.68761335   8.97642389   9.18997755\n",
        "   9.31866582   9.36587056   9.34740836   9.28893189   9.22171529\n",
        "   9.17751587   9.1833565    9.25708583   9.40444579   9.61812821\n",
        "   9.87897556  10.15912843  10.42660281  10.65054491  10.8063004\n",
        "  10.87946503  10.86825119  10.78378163  10.64826203  10.49133265\n",
        "  10.34519853  10.23933827  10.19566084  10.22490593  10.32487947\n",
        "  10.48081414  10.66779556  10.85485568  11.01006072  11.10575781]\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Draw a plot to compare the true relationship to OLS predictions. Confidence intervals around the predictions are built using the ``wls_prediction_std`` command."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure()\n",
      "plt.plot(x, y, 'o', x, y_true, 'b-')\n",
      "prstd, iv_l, iv_u = wls_prediction_std(res)\n",
      "plt.plot(x, res.fittedvalues, 'r--.')\n",
      "plt.plot(x, iv_u, 'r--')\n",
      "plt.plot(x, iv_l, 'r--')\n",
      "plt.title('blue: true,   red: OLS')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "<matplotlib.text.Text at 0x469b110>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEICAYAAACpqsStAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYlFX7B/AvIiKoCIYiCYqiKKKAC2qaQsninqmYS2pu\nmOXSpr8SFypNM99S0XrTNDH3pUxBETdw316XyA0FV1BTGBBkZ+7fHydG0WF75pkN7s91zZXOPHOe\nw0g3h3Pucx8TIiIwxhgzSlX03QHGGGPScRBnjDEjxkGcMcaMGAdxxhgzYhzEGWPMiHEQZ4wxI8ZB\nvJJzcnLCgQMH1L4WHR0NR0dHHfeoYgoJCcGIESP03Q1WAXEQr+RMTExgYmKi726gSpUqSEhI0Hc3\ntKa8n3FqaiomTpwIe3t71KhRA+7u7lizZk2Ra0r6AfzNN9+gSZMmqFWrFhwdHTFkyBCpXWcGrqq+\nO8BYoZL2neXn56NqVcP4dtV2X3Jzc+Hr64v69evj5MmTcHBwwP79+zFq1CgoFAp8/PHHAIr/ARwW\nFoZ169bhwIEDaNy4MR4+fIhdu3Zprb9Mv3gkznD69Gm4ubmhTp06GDNmDHJyctRe9+Jo+b333sOs\nWbNUfw8PD4enpydsbGzQpUsXxMbGlun+3bp1AwB4eHigVq1a2Lp1K6Kjo+Hg4ICFCxfC3t4eY8aM\nQVhYGLp27Vpsn3JycvDZZ5+hUaNGqF+/PiZOnIjs7OxyfRbFcXJywsKFC+Hu7o5atWpBqVTi5MmT\n6Ny5M2xsbODp6YmYmBjV9Tdv3oS3tzesrKzg7++Px48fl/lev/32G+7evYutW7eiUaNGMDU1RUBA\nAJYuXYrZs2cjIyOjxPefPXsWAQEBaNy4MQDAzs4O48aNk/aFM4PHQbySIyJs2LABUVFRiI+PR1xc\nHObOnVum9z4/Ejx//jzGjh2LlStXIiUlBRMmTEC/fv2Ql5cHAPjwww/x4Ycfqm3n8OHDAIC//voL\n6enpCAwMBAA8fPgQCoUCd+7cwYoVK0ocqQPA559/jhs3buDixYu4ceMGEhMT8dVXX5XpaymLTZs2\nYc+ePUhNTcX9+/fRp08fzJ49GwqFAosWLcLAgQORnJwMABg2bBi8vLyQnJyMWbNmISwsrMio2cPD\nA5s2bVJ7n3379qFXr16wsLAo8vyAAQOQnZ2NEydOlNjPTp06Ye3atVi0aBHOnj2LgoICDb9yZsg4\niFdyJiYmmDRpEho0aAAbGxsEBwdj48aN5W5nxYoVmDBhAry8vGBiYoKRI0fC3NxcFXCWL1+O5cuX\nl6vNKlWq4Msvv4SZmRmqV69e4rVEhJUrV+L777+HtbU1atasiS+++KLYQFleJiYmmDJlCho0aABz\nc3OsW7cOvXr1Qo8ePQAAvr6+aN++PSIiInDnzh2cPXsWX3/9NczMzNC1a1f07du3yA+hixcvFjtP\nnZycDHt7+5eer1q1KmxtbUsd1Q8fPhyhoaHYu3cvfHx8YGdnh4ULF2rw1TNDZhiTjEyvns9Aadiw\nIZKSksrdxu3bt7F27VqEhoaqnsvLy8P9+/cl96tu3bqoVq1ama599OgRMjMz0a5dO9VzRASlUin5\n/i96/nO6ffs2tm7dWmSuOT8/H2+++SaSkpJgY2NTZCTdqFEj3L17t0z3sbW1VftvkJ+fj8ePH8PW\n1rbUNoYNG4Zhw4ahoKAAf/zxB4YPHw5PT0/4+/uXqQ/MePBInOHOnTtF/vzqq6+qvc7S0hKZmZmq\nvz8foBs2bIjg4GAoFArVIyMjA++8847kfr24aFejRo0i93/w4IHqz7a2trCwsMDly5dV909NTcWT\nJ08k37+k/jRs2BAjRowo8vWmp6dj+vTpsLe3h0KhKNLX27dvlzlDxdfXF3v27CnyfgDYvn07zM3N\n0alTpzL32dTUFIMGDYK7uzsuXbpU5vcx48FBvJIjIixfvhyJiYlISUnBvHnziv0139PTE+vXr0dB\nQQEiIyNVc9kAMH78ePz3v//F6dOnQUR4+vQpIiIiSl2EK2RnZ4f4+PgSr/Hw8MClS5dw8eJFZGdn\nIyQkRPValSpVMH78eHz00Ud49OgRACAxMRFRUVFFrnm+z5p49913sWvXLkRFRaGgoADZ2dmIjo5G\nYmIiGjVqhPbt22POnDnIy8vD0aNHER4eXua2R4wYAQcHBwQGBuL27dvIy8vD3r17MXXqVHz55Zeo\nVauW6trc3FxkZ2erHvn5+QgLC8Pu3buRnp4OpVKJPXv24NKlS+jYsaMsXzszMMQqNScnJ1qwYAG1\nbNmSrK2t6b333qOsrCwiIjp06BA5Ojqqrj179iy5ublRrVq1aMSIETRs2DCaNWuW6vXIyEjy8vIi\na2trsre3p8GDB1NGRgYREb3//vv0/vvvF9uP//73v2Rvb0/W1ta0detWio6OLnLvQvPmzSNbW1tq\n2LAhrVu3jqpUqULx8fFERJSdnU0zZsygJk2akJWVFbm6ulJoaCgREd25c4esrKwoJSVF8ud04MCB\nIs+dOnWKvL29qU6dOlS3bl3q06cP3blzh4iIEhISqGvXrlSzZk3y8/OjyZMn04gRI1TvdXNzow0b\nNhR7v5SUFJowYQLZ2dmRhYUFtWrVilatWvVSn0xMTIo8Zs2aRb///jt16dKFbGxsyMrKitzd3Sks\nLEzS180MnwkRHwrBKr7169fj8uXLmDdvnr67wpisOIgzxpgR4zlxxhgzYhzEGWPMiGklT9wQCiox\nxpgxKu8Mt9ZG4kTED5kec+bM0XsfKsqDP0v+PA35IQVPpzDGmBHjIM4YY0aMg7gR8PHx0XcXKgz+\nLOXFn6f+aSVP3MTERPL8DmOMVVZSYiePxBljzIhxEGeMMSPGQZwxxowYB3HGGDNiHMQZY8yIcRBn\njDEjxkGcMcaMGAdxxhgzYhzEGWPMiHEQZ4wxI8ZBnDHGjBgHccYYM2IcxBljzIiVGMTHjBkDOzs7\ntG7dWvXctGnT4OrqCg8PDwwYMABpaWla7yRjjFUo2dnAtm1ASorGTZUYxEePHo3IyMgiz/n7++PS\npUu4ePEiXFxcMH/+fI07wRhjlcKTJ8CMGUCDBsCPPwL//KNxkyUG8a5du8LGxqbIc35+fqhSRbyt\nY8eOuHfvnsadYIyxCi0/H/jpJ8DFBbh/H/jf/4CDB4EWLTRuWqPT7levXo2hQ4eqfS0kJET1Zx8f\nHz4BhDEdiYg4jKVLo5CTUxXm5vmYMsUfvXt303e3Kre//hLTJ3v2AG3aqJ6Ojo5GdHS0Rk2XerLP\nrVu30LdvX8TGxhZ5ft68eTh37hy2b9/+cqN8sg9jehERcRhTp+5FfPw81XPOzsFYsiSgbIFcqQTS\n0oBatYCqGo3xmAQ6O9lnzZo12L17N9avXy/l7YwxLVm6NKpIAAeA+Ph5CA3dV/yb8vOB778H+vcH\n6tYFGjUCzM0BOzv11xMBOTky9pppotxBPDIyEt999x3+/PNPVK9eXRt9YoxJlJOjfvScnW1a/JtM\nTcUC25Ah4tf+J09EYE9IUH/99esiwA8eDKxbJ0uGRYUREQHMmaPTW5YYxIcOHYrOnTvj2rVrcHR0\nxOrVqzF58mRkZGTAz88Pbdq0wQcffKCrvjLGSmFunq/2+erVC8Q0yaNHL79oYgIsWCCCeIMGz56r\nUUP9TVxcgGvXgB49xDyvkxPwxhvAjh3yfBHG6NEjYNgwYMoUoJtu1x/4tHvGKhB1c+Kujadj25uZ\naLlrK/DNN8DYsfLeNCsLiIoSgb9fP3nbNnRE4reRzz4DRo4EvvwSsLSU3JyU2MlBnLEKJiLiMEJD\n9yEvywR9ks/j/X9Ow6Lza8DcuUCrVvruXsWyeDGwZg2wahXQrp3GzXEQZ4wJubmAlxdgYwPMnw+8\n9pp++pGXB8yeDbz/vlgwrWgyMsQisJmZLM1JiZ2cQ8RYRVStGrBxI+DqKqY59CU3V9y/bVsx5164\nW7GiqFmzzJdqK3+fR+KMMe179AhYuBBYvVrMHX/+efEpjIYoOVk8XFwkvb2s+fs6yxNnjLFyqVsX\n+O474NIlsaHoxAl996hsCgrEnLebm0gflEhS/n4Z8XQKY8Zu0yagdm2gZ09996R09esDS5bouxel\nUyqBrVuBr74Sn214ONC+veTmJOXvlxEHccaMVX4+8MknwO7dIuAYuydPgBs3xPy5vgUEAOnpYier\nv7/G6wol5u9riKdTGDNGmZnAwIHAlSuiIt5zRZWMVlwc8NZbQPfuolCUPtfVVq0SUz4BAbIsDE+Z\n4g9n5+Aizzk7z8DkyX4at80Lm4wZm8ePgb59gaZNRbCpVk3fPZJPXh6wZYuYP8/OBgYNEpuTGjeW\n/15EYrHS1lb+ttUozN/PzjZF9eoFmDzZ76XsFM4TZ6wyOHlSzNF+/bV+0we1iUh8nTt2AAMGAB07\nytfumTPA9u3A778DTZoAe/fK07YGCtMPo6LmcRBnjFUSISEiT9vREWjYUPzX3l4U9HqRUgmMGSMO\nYrC0FFNRAweKaSg9/yAsmn7Im30YYxIY5UEStrai0uKpU8CdO8Ddu2KqSaF4uXhXlSqiSNf06UDL\nlvrpbzHUpR+WBwdxxio5dRtR4uPFIpxBB/JJk15+Lje3+MMsRo3Sbn8kKi79sKw4O4UxQ7Z4sZi/\n1SJtbkTRuWrVxKjbiBSmH/6MIEnvN66vlrHKggj44gvg559FISst0uZGFFa6wvRDF8RJej9PpzBm\naPLzgQkTgL//Bo4c0XoKnDY3orAXBAWJfHhLS2D9esDGRjVlZfneZuBx+ZvkkThjhiQrS2RNJCYC\nBw7oJIdZmxtR2Avi4oCYGLGZyclJtaGpd+9u6HD9rKQmeSTOmCFJSBDV/ZYt09kmnsKRYGjorOc2\novQw7EVNY0Mk0hsvXhR/b9RI5Kc/l95YUMtaUtOcJ84YAwAoxwch60IcnpIlEuZugFVDa9SpA5w8\neQTVJk3BqxkK5FatCsXyJQh4p7e+u2u4np8y2bABsLYWZ29GRopj3PbtA1auBNW2Rlyc+IVr/37g\n0CEgNZV3bDLGyujRI5FiffIkcOqEEt9Gd0RbpfiV/lb1FvjQ/nccVzgjNbUKDsEXPogBAOyo1gVP\nVi7CyJGdSmxf7txzXeeyS76fj4+YMgGAwEBRRiAxEbCzA5lWxd69wObNInADgK+veLz5JvDqqxJi\nJ2mBlppljGkob/R4uuvsTYcse1JDKwX5+RHtfXMB5dS0IaWFJRFA5OBANHMm0dWr5O8fTABRBHoS\nAXQKXrQDfWmLyQCa13glLf/8DikGjyfy9ibq2ZNIoSAiovDwGHJ2nkFiHkE8nJ1nUHh4TIn9Cw+P\nIX//YPL2nkP+/sGq66W2J5VG9+spPivy8lJ9Hnl5RBs2EHl4ELVqRRQaSnT1KpFSWfStUmInB3HG\n9EGpJPryS6K1a4u9pLiAJkVmJtH48TcoxuR1VVS618VHvHj9OlFSkgg4gYGqwENE5O09hwCi2lDQ\nJgRSbSjIHok0z6U/3fcZQunVX6GnsFS1mTcgkIhIFfxffAQEzCzx6y0ucEppTxOS7pebS/TDD0Tx\n8arPMTOT6McfiRo3Jnr9daLw8JcD9/OkxE5e2GRM17KygNGjgVu3gD/+UHuJXLsos7JEqvnyLxPR\nKWsHMqgWAOA0vDAhsS3mRhwu2t6WLUXeX5h+mAZrDMEW1Z8PN26FGZFfA0ol6PWuwInjuGHdHr2P\nrkDQf4DMTAu1/Skp97z4TUezdJ7LXu77HT0KTJwozg8dNAjpq7Zg2TJx/kWHDsBvvwFdumilq5xi\nyJhOJSUB3bqJIk3R0aJgkxqa7qJUjg/CvWY+OGPtB6dln+Fshgs65tzEMKzHZgTCH1G4cOs/pbZX\navphlSow2R0BBAai6c192LrPGqdOAadPTwUgdiEegg8i0Au1kVpi7nlJgVPXueyl3i8oSMx9+/oC\nw4cDQ4cCc+ZAGbEHa/Y7oHlzkea/fz+wc6f2AjjAKYaM6c6ZM6Ks6sSJYjdmCdXzNBl5JiUBj7Ze\ng0faYTgAgDIeQ9qNx+ZTiwFANaIuS3tlSj+0tlaN4N3//eOPP17F//2fGXwyouGC6wCA9TVeAyb/\nXOy9Sgqckyf7Iz4++IWDhmdg8uQeJfZfqilTSrlfYb43IA5PvnwZx2NrYWonUbplxw4xAtcFDuKM\n6UrNmsBPPwF9+pR6qdSRZ0SEOEPhkmmSeKJ1a+DwYSjeWSSpPUAE8vJmgXzwQXs0anQY6cMIeALk\nwRS1m7RB57bNi31PSYFT17nsJd0vIuIw6v59Dx0A3KhZD6c/WoaI92shJgb49ltg2DAdV7ct9yx6\nGWipWVZJybnAZyzUL/J9UezXnp1N9PHHRA0bEh0+TEQPHhRZpCxve7L5d7H04uYrtNl+CilM69Dt\nodPFIqAa4eExFBAwk7y951BAwEyD+7cu/ByfX+itUuUpvfPOLUpP17x9KbGT88SZQVO3wOfsHIwl\nSwIq/I7CshznBYizhYcMARwcgNWrgTp1NGtPW4iAXT/dQ2zwJvzP51P8aj4BtR+8sCnG0BRu3MnL\nAyIiEPDOIkRFzX3psoCAWYiM/FrjXHZJsVPznx0v01KzrBLSdWqZLE6fJvr8c+3fZ/x4eujqTfvM\netCKhYoSU9cMSVYW0TffEB2t6v3sHzQwUN/dUq9t22d9HDiQ2rX7Ue33o7f3HFly2aXETs5OYYYl\nNxe4fVuczkJkXGVSL14UE6J9+gAeHlo9rZ0IuLUvDvWuxMA3LxLjzwTp+5SxMqteXazrtu9mCQA4\nV9ULP7iuQM62XWJnoyE4fx7o1w+4dAkAkOfZHh/X+gWxsSPUXl69eoHe6rJzEGeG4fPPRVGgWrVE\nPlbjxoC5OV5LV19j2dY0Q6tBslyOHwd69wZ69gQ8PcWv30OGaG11Kz8fmDCuAGaJt8QTrVsDK1Zo\n5V7aZL59AxAYCKsTUThxxRo/jT6NbJfWULZwBTp1Anr1AlJTtd+RwnTBwvutWSP+PX19kXX1FuI8\nA+F6dx+otjXWrv2r2JRLfQ04eE6c6U5uLpCSAtSv//Jr164BZmbisFszM/FcTg527zmKKZ8dfClj\n4XjtA6h3Ox5o314cmlD4KCbvWqtCQ/FX3B3MuGKKjPzqWq3rkZEBjBiYic8uvouOzo9QtZ4t8Ouv\nhjmfXAJ1c8d2dt2w4LPHWHKsPRrk3xYXdu0KHD6s3c68WOtk5Uo8yTLDf9da4ocfRKrg/PnPjuYs\nbm0hIGBmifPlZcFz4sww3blDFBxMZGdHNGdOud9ebMZCYiLRjh2izkdAAFGdOuJe2pKdXWz/dFHX\n4/59os6eT+lGvU5UMOxdopwcWdvXldI+r0deovZIollDuuD7CT158tybx79cp0USpZLo77+JCgqK\n1Dp5fENBs2YR2doSDRtG9Ndfmn5d5csAkhI7OYgz7Tl1imjAACIbG6LJk4muXNHu/ZRK9YUp8vOJ\nxowhmjWLaNUqogMHiG7cKD4IKpUiQNy+TbR9O9GUKUTu7kQdO6q9XBeLr1evivobX4YoSblzV8kF\nOHRISvpnqZ+XQkHKwEDat1VBb79NZG0tYveZMyQCeOEb+vRRfQ7h4TEU0cCTztduRKdecabITeFF\nb1oY/Lt0IfrkEyIXF5GPeecOkUJBT/sE0hcTFWRjIy69fl3656FJiiQHcWY4MjOJWrcmWr6cig6l\n9CAri2j1aqLZs4lGjhT/MzdqRFS3rvrrHz0isrIisrcXo7QFC4hOnCg2t7mwSJS6jAWNjR9PqW28\naX+1nrR+uQYjTy2Q+htIeT+vpCSiefOInJyIjtQSo+aCVx3Ev5+lJaU5NaE9NZrTbTiqGguv0aJo\nPxo0eHaj5s2p4ORpOnlCSUOH3iIrq7tUtWomNWp0lNasOSHjJ1R+HMSZYSko0HcPSpafL0sz2hyJ\nP27tbbBpeFK/bqnvKyggOrBdQccaBJK9hYJcXIhGDXhC/RyX0SBsoatwUZXLrQ2Fqr3cXKKcrm8Q\nAfS4iRdNeEdBdesSNWyYQdbWMVqfBisP2YP46NGjqV69etSqVSvVc8nJyeTr60vNmjUjPz8/UqiZ\nl+IgzioTbe2GXL+e6JCZ30u1qQ2F1N9A5Pi8cnOJYmNFJV8Hh+MEFC2XCxCZmaWThQWRqSlRw9oK\n2mkRSIP9FfTjj0S3bhnmHgQpsbPE2imjR4/G5MmTMXLkSNVzCxYsgJ+fH6ZPn45vv/0WCxYswIIF\nC8q3msoqlvR0sevO1ABzt3VAG3U9li8jZAbPRT+PHKBxoEghNLAMFKn1XeT4vMzMgFatxGPdugjc\nu/dakXK5ANCly3JERPwfLCwAExNrAFvQ97k2jGoPQklKi/I3b94sMhJv3rw5PXjwgIiI7t+/T82b\nN5flpwkzUufOETVpQhQRoe+eVAhKJdFXIQUUZjWJsl09REqKGoZQT0Zv9Vhk6kelGImr8/DhQ9jZ\n2QEA7Ozs8PDhQ7XXhYSEqP7s4+MDHx8fCT9imEE7dUrsalu6VGyUYBpRKoFpU3Ph+9sodG+ZhGqR\nMUDt2i9dJ9eBEZrSdWVBuftRarlZHYiOjkZ0dLRGbZS62efWrVvo27cvYmNjAQA2NjZQKBSq1+vU\nqYOUlJSijfJmn4rv6FFRG3vNmkoVwLVyWG9QEJTX4hB7ozrynubCo4sVzLZtBCzUn44jx6YSJuil\nKJhCAezZA7i6Am3aFHlJSuws90jczs4ODx48QP369XH//n3Uq1evvE0wY3fihAjg69cDfn767o3O\naGsEnH8lDlWPxsADQMFrXWD65zZxskAxKsxcrgGQUitdsps3gcWLxVltXbsCn34qS7Plrp3Sr18/\nhIWFAQDCwsLQv39/WTrCjEjTpuJsSAkBPCLiMAICZsLHJwQBATMREaHlLdUy9kMbBY4SE4GTF0Uh\nKGV7L5juDi8xgAPSFxTLhEgsVOfkaN4WE5KSRC2d9u1F9a/YWODPP8UxfXIoacJ8yJAhZG9vT2Zm\nZuTg4ECrV6+m5ORk6t69O6cYsnLT1fZ0bfVD7k09Fy4QOToS/TBH7FAsawqhrAuK+flicXrJEqJB\ng4jq1yeytCT68EP11//5J9G4cUS//abdEgcy0vsi8JMnRN9/T5SWVuqlUmInb/ZhOmMo2QC63qSi\nzu7dRE1fSaHNm6V9DbKdgBMVRdSihdhrvnYtUUJCyVv6b94UAX/AAFFgpHFjovfeIzp5Utr9tcxQ\nBg5lJSV28hmbTGcMZS5Xaj/kymb4+SclsqfNwv+cD8Mq8DCA8peslW0u188PuHKl7Nc7OQFTpoiH\nUineW1gB0AAVPwU2S/658AcPxPxYu3bytlsKDuKsZGfOAJs3A4vUH7RbHlqdy9VBPzRKqQsKAsXF\nIe62ORo/skDX5o9gEblDNyfqJiSIE3xnzhSlfuVSpQrg5iYexSko0OsmMJ0MHHJygCVLgIULgeBg\nnQdxnk5hxUtIEEWgduyQpTlj3xyiiezXvFU3y2/gKIpyadvVq6Lg1yuviHK9KSnav+fzHj8Whadm\nzBAHN+uBVqfwlEqxRuDsLCoqxsVp3KSU2MkjcaaeQiHyv7/4AnjrrZdelpIvbeybQ6Tatg2wPWsG\nHwDUoAFMYy+KLAVtyc8XI+/Fi8W0x40b+tmy/8orwIEDYjOYqyswfjwwfbp4Xke0uqFn4kQgOhpY\ntgzoobsNQi/R+EeHGlpqlulKdjZRt26i7rIaxrZYpC9Pnog1v6ZNic7sU4jRmi6KWMXHi3vdvav9\ne5XVnTtEEyaIgztk+s2urGRbBH7RtWuyH8whJXby8WzsZbNmiQWrLVvEvOcLeMdg6U6cAN59F3jj\nDTEgrllT3z0yEAkJgLk50KCBvntikHSyY5NVAtOniw0nagI4YDhZJqXRyhb5kgQFoeBqHBKSLDAy\nbSO++9kaAwZo73ZGqUkTffeg/BISxNmtxZRB0DcO4uxltWqV+LLcWSYKhRj4X7kCtFwShFpJccgz\ns8TaHhtQrZ41rK2BtyKC8EpKHExrWgIbNsC2aclzvLouEpWRASgOxMExIQbNAMR6v4PqA/aW+j6N\nf9BcvAi4u+smy0WbLl4EoqKADz4AatTQd2+EBw/ECcnr1okdlq+/ru8eqSfrhM6/tNQsMxCaZnek\nDRlP8Q296WSdnuRST0GNLR/QdOdttLXNPMqoaadq9B/HNjR/PtH//R/RNXtv1fN/Vh1AjR3z6O23\niebOJbrTYzzldPYucniurjYWpaYSzf2qgEZZ/U4ZVa3ETZo0KZIJUtyOQY3XFkJDRfZQUpKsX5Ne\nXL9ONHiwOEx74UKijAz99eXxY6Lp08X8/Ucf6TSzRkrs5CDOJCnvYpFSSXTkiNjod7SqtypqPe0d\nSMo9kUT9+on/cVq3Fq+5uIg3FHruRPKCRf+hAgtLetykPZ10H093LZqp2jvuGEjLlhG1a/eTrFvk\nixg/nnK7eNM1557Uw+oo3a3tRplu7cRW9EGDiixelhSoJf+gyc8XB0+7uoo00IokNlYcQ2dnJ842\nLcNWdVnduiWC9/vv62VhmIM4K7/kZJGFItN5ky/KyRGxrZtnGk20/52WLSPK83sWkF/K1lAoxP/E\npT2fnk509CjRsmVEDg5EAOXYN6TfQhU0ZgyRpeU/BBD9jPF0CN4UgZ5UGwry958l6etISyPas0ek\nPJ+r7a2KuOl+/cULxWxVLylQS6rF8uQJUe/eRN27G9xxbbKKjRWpPcnJur/3vXu6v+e/pMROnhOv\nzPLygMGDxZyqzLvqaHwQkmLicPeWEjZ1m2Pfk+0w6+kHk6A+wPANQFCQ+iPHrK1FVsyLXny+Zk2g\nSxfxGD4cCApCtRUr8K61Nd4FMGDAFUyatAIut+LgA7EtPArdMe3gYni3fAQ7N1tMi5+AhtlxgKUl\nzn6yATnVrVFQIDYZtlsRhJpJcVBkW+Ab8y9hc+sC4jq8iw4+lnBwsQTOAPDyQs0tv5aYg13SIrCk\ntYX33xeLbD/+KM4oq6hatQJ+/VX9a0SarwFcvQpYWQGvvvrya8aWOaOFHyY8EjcWkyYR9egh+yj8\nn3+ILlkNvLPHAAAfKElEQVR3fja0dHXVy+gmPDyGTtk6EwH0yLwmJXbxpvx2XpRXy5pyLGtTYv02\nqj4etg+k/v2JBg4k+r3FF/Sk2ivPdlhWt6T84SOezT0X99uCGiWNxCWtLaSmllygqjIIDydycyMa\nNozom2+Idu4UufEFBeqvz8gQu1cPHyYKCRHvffVV0Y6BkRI7OYhXVj/9JKrXpabK2uz+/WKn9cM6\nzUVUatdOv7/2qwu4SiXRo0dE/v7qp3UOHiRq21a85u6uUf9LC9Ra24hSkeXlEZ09S7RmDdG0aWK9\nxNGRKChI/fUrVhA1a0bUubNYqDx2rPiAr2dSYidv9qmMTpwA+vcHjh0TBzzIIC8PmDMHCAsTD9/2\nqcVPmRiK1BL6WNJr5aSXI8AqIz0X25KDlNjJQbwyysoC4uIADw/N2gkKAuLikAlLvPV0A8zqWmPN\nGoBP7JPJnj1AQECxm65YxcNBnOmWj4+qlvQ1j0A0O7eF440ciIDPPwd27QKOHNFpwSimX7ztnunU\nw1Rz2AHIq/0KmkevkHBiK3tJfj4wYQJw6RIHcFYmHMSZJH+GpcLhUiqs7Z1g/tcZw533NibZ2cDQ\noUBmJrB/P1fNYmXCY6fKYPNmcYK5THas+AdNxr2BhoM7wfxePGBrK1vb5SXl1HqDNW2aqDO+axcH\ncFZmPBKv6DZsAGbMALp1K7WwVVlsWpsLzw/fhM3YAaj705d6Lbyk6yJXWjd3rvg34oUFVg68sFmR\nHTkCDBwIHDwodsBpaO1acdDP4R//hvNbmrenKa5rzioaXthkz8TFAYGBooymJgH83zTCuymWWPB4\nAw4ctIZzC/0HcEB7dc11XoecMQ1wEK+IkpOB3r3Fr+f+/pq1FRcHxMTAEcDJHkGwaqGmromeyF3X\nHNDhFM3evYCvr9FvTmH6x5NvFZGVlTgod9w4jZtKyRanmWS6ecFq4wqN25PTlCn+cHYOLvKcOATX\nT3KbS5dGFQngABAfPw+hofskt1kEETB7tjj84MEDedpklRqPxCsiMzPIcS5YwplkpJ/9B/kd+qDe\n3t8MLo1QG6fWa/Xoudxc8YP12jVR+oC3tjIZcBCv5Iqb/1UkZiKtW1/gTV/Ui/pWlja1oXfvbrK2\nrY0pGgBiWmroUMDJCTh0CLC01Kw9xv7FQbwSK27+Nz+7ALYTFsPcqQnaR86XpU3AONL+pkzxR3x8\ncJH+iymaHpo1PGcOMHYsMHGi8Z+HyQwKpxgau7w8ICREnFBfu3a53qo+RY8QZhkI9xppaH07AqYW\n1WRo07jS/rRSdVCOgwxYhccphpVNfr441SYrC7CwKPfb1c3/uuESXLMT0OxGNCIPniz3tIhW55R1\nRO4pGgAcwJnWcBA3VgUFwHvvAWlpwJ9/AtXKN2IGns3//owguECUlB2GDfi/1wfj03MXJE2LaG1O\n2Vjk5IjaJzY2+u4JqyQ4xdAY5eQAI0cCSUnAjh2i3oYEhSl6LhDnUPbCHvxm8QY+nd5ZcqqdJml/\nURt2YprXu/jepQ8WuPdXXweloEAUijI02dni3MtmzYDQUH33hlUiPBI3RsuWiXS18HBJ0yiFCkfU\n1UduA1KAvy2boNrquQjo3Q3ffXdQ7XtKmxYpd9pfSgrwzjvIPnsOrz15Cktle1xFC9yHPaZO3Vuk\nTQBAfLw4zKJxY6B162cPDw+R+aFrWVnAL7+IvPw2bYBt24AOHXTfD1Zp8cKmMcrPF0WSZCiUlJdL\nePvNNMxPDkLrE8+OItPZAmVBARAZiXe/jcT6I0sBFJ07Vnu/nByRax0b++xhbQ2sXy9fv8oiPx9o\n2VI8Zs0C2rXT7f1ZhcMLm5VFVfX/bOXOz87Px9WWg+BUbzZa/r0FeG6QLWuq3YMHwMqVopZLixZF\nXzM1BXr3xr3vzuDFAA4UM/I3Nwfc3cWjNHv3Atu3A+3bA25uIuDKNV9dtaooMmZnJ097jEkgOYjP\nnz8f69atQ5UqVdC6dWv8+uuvMDc3l7NvrByk5Gdf6/UxFEnZ+PqE+0slPGTZDXnxIvDdd0BEhAjg\nJXx/aG1BtHlzUQDs5Elg1Srg8mVRq3vmTODDD1++/tIl4OpV8WelEvjrL+DoUbGIPGrUy9dzAGf6\nRhLcvHmTGjduTNnZ2URENHjwYFqzZo3qdYnNshfl5hLNnEn06FGpl/r7B5NIRi76CAiYqfb62zN+\nomumLejvY6ly95ooJ4do9myiunWJFi0iSkkp9S3h4THk7DyjSN+dnb+g8PAYefumVBLduUN09676\n19evJxowgOjtt8UjOJhozx6itDR5+8GYGlJip6SRuJWVFczMzJCZmQlTU1NkZmaiQYMG8v50qezO\nnBFnLTo4lDiCLVSe/Oy0Pw/BYsEcxH5/DL07l2+DUJkkJ4tt5hcuAK++Wqa3aKMOilomJoCjY/Gv\nDxsmHowZCUlBvE6dOvj000/RsGFDWFhYICAgAL6+vkWuCQkJUf3Zx8cHPj4+mvSz8njyBAgOBrZu\nFVMR775bpo0iZZqOCAoCXb0G82Nn8Hu/TRg9talcvS7K3h7YuLHcb9PKJhvGDFh0dDSio6M1a0TK\nkP/GjRvk6upKjx8/pry8POrfvz+tW7dOo18JGBGlphI5OhKNHUv0+HG53lradER4eAxdsGmoevFe\nFx9tfAWMMQ1IiZ2SRuJnz55F586d8corrwAABgwYgOPHj2P48OGa/USp7GrXFkepNS3/CLmk6YjC\nRc+lCjd44A5OwwsTEttibsRhzUe+Bw4APj58uAFj+iLlp8WFCxfIzc2NMjMzSalU0siRI2nZsmUa\n/TRh2lO46FkbCtqEQKoNRYmLnmWiVBJ9+y1Ro0ZE9+7J1lfGKjMpsVPSbhEPDw+MHDkS7du3h/u/\nubpBQUEy/mipwHJzxXz3vHmlXyuT9HQrAEAarDEEW5AGsaFHclEqpRKYNk2cnHzsGMCL2ozpjeQ8\n8enTp2P69Oly9qViu3MHWLFC5Cq3aAFMmqST2+bvCEeLvyxwQs1rknKw8/JEXez4eODwYaBOHY37\nyBiTjgtg6ULXroCnJ5CeLua8Dx0CBg7U/n0vXEDW0NFQOjZHkyYzi7wk+SzKOXNEvZN9+ziAM2YA\nuHaKJjIygMTEZ48331Q/tRAfDzRqVOx2ea24fx8ZrToi2HwRvroyGEePynTQQVqaOFrMzEz+PjNW\nyUmJnRzEX5SWBly/Dty/D7Rtqz4oBwUBmzeLqYUGDZ49vvhCbPHWp6Ag4OpV5J+7iKXKSeh1bt5L\n5UoYY4aJg7hUS5eKoHz9uigt2rSpCMozZwKdOr18fVKSKAFrbW14J7b4+AAxMQCApC6BePXoFv32\nhzFWZlzFsDQ5Oeq3sHfqJOasmzUD6tcvPTCXcSu5PuSbW6IqgKRX2+HV8BWaNZadLaaAdDkNxBgr\nl4q9sEkk6nfMni2mOYrLpunQAejWTWwXN7SRdTnk5QHDaANOOwXC/u/9qtrgkuTnA++8Iw6gYIwZ\nrIo7nfL338DHH4spksBAYMAAoGNHWQ5SMEREwJgxwD//iBPbNFp3JALGjROLtTt3Sjq/kzFWfjyd\nUigzE3j7bWDKFOD99ytFJsXMmaJU9sGDMny5M2aIH4IHDnAAZ8zAVcwgbmkJXLlSIedyXzy95zuX\nHPyT2grbTo/CsWNAjRoa3uCHH8RQ/sgRcXgCY8ygVbwoV6iCBvCpU/dievwjuCAOtZGKuvtuY7TV\ncUSeB2xtNbxBQQFw7pw40kzjxhhjumDckU6pBNatE0X8K1DQLu6szKVLoxAfPw8u8IEPRBphJAXA\nvOVBNG7sqvmNTU2B337TvB3GmM4Yb+TLzgZGjwZu3wb69dMsE8OAlHRWZuHpPbZ4BAD4Gy0xBJvg\nWW2x7jvKGDMIxpmqkZwM+PmJX/8PHKgwARyAarT9vPj4eQgN3Qdz83zUQTLSUBvh6IXXcQxpsNb8\nMGHGmNEyviB+4wbQubN4bNokdk5WICWdlTllij/M7XfidRxDX0QgDdbSC1kBIhc8N1eD3jLG9M34\nplNmzwY++UQcIlwBlXRW5pMn3ZCV1QkdOqyEhUWS5ocJf/KJSGeZP1+DHjPG9Mn4NvsolRV2ww7w\ncgZKJiwR7NQCHfw/wO7dTbFnj0w1tn7+WaQTnjxZoaajGDNmlWOzTwUO4MCzszIdR4yEu+I2AMCm\noArei26Ko0dFRVuNHTwofqM5epQDOGNGzviCeCXQ29MZqJYNALhfuwWCbdfh6F6gbl0ZGr9+HRg6\nFNi4URT8YowZNcMe1t69K1IJDVxExGEEBMyEj08IAgJmIiLisLSGCgpEwSlPT6T2G4kY24H4oM0J\n/BljLU8AB0T7X30lDrBgjBk9wx2J//MP8MYbwH/+A7z1lr57U6yS8rpLXXAMCgLi4kSZgA0bgLlz\nQadPY9PEw5jykys+/xzY9pHYgyObH36o8FNSjFUmhrmw+fSpCOA9eiCio6/a3YuGIiBgJqKi5qp5\nfhYiI78u+c3PHeCAwEDcnrkSYz+qhYzMKlizBnwiD2OVTMVY2MzPB4YMAVq2RESH7pJGucVtW9eG\nkvK6i+0HkahbbmkJACAvL6x5bQWmd6+NadNE5l8FqiLAGNMiwwoVRMCkSeIEnpUrsbTPl8XsXpxV\nbFDWaHpDguLyutPTE19KFfzPlWQ08d0I1+OHgIgIYMMGpL4ThNG5K5C40RoxMUDLlrJ3kTFWgRlW\nEM/PB2rXBrZtA8zMShzlFqf4besi8Ms9Sp8yxR/x8cEv5XUTVUN8/Dx4ogM64AwAwPduVUTvbQmX\nTasxZ3UWfvwpBU+f/gInp2P47rsaaNlS5h8yN26IYf0ff8g8sc4YMxSGFcTNzIBvv1X9taTdi8Up\nbXqj1FH6i4uNz+dRq3mtd+9uqJKXC8fRY9Aq9S4AwNMhC8NMvQEAlngKAEiAE3wQA7PqV5ExqC3S\n0h4iJ0ckfcfF9cQnnwTD1FTG3xZSU4G+fYGpUzmAM1aBGXSawpQp/nB2Di7yXGm1QkoK/CUVl1KJ\nixOLjXv2iKD9vOdfa9YMaN4csLZGz8G90Mr23wOYvbzw6q4/VP14HcewGYFoi/O4i4bIyqqHxo1/\nVQXwYvuhicJ1BV9fcbIRY6zCMqyR+AsKR6WhobOQnW1aplohhdMbzwdrEfh74LvvDgIAfkaQaupj\nKhbjtTtXxA7GS5eAs2fFm+rXB1a8cFr8vwuRcHYGFi8W/7WzA2xsgLQ0UFAQ7s5cgZhd1qhS5UOY\nmT1EWp4dhmALAKBRo2+wfPnr+O67h2r7XtI0UZkRiSkUpVKkEzLGKjT9B/GCghJ/3e/du1u5phhK\nCvxLl0YBAFwQpzpUwQKZMEt5CKAlMHgw8OmnwPffA7/88vKW9A0bkDc2CAnTV+Cmwhq3DwO3bomS\n5rdvWyMhYQvIH+jaFejZ0x49epxDRMSPyM01+bcffkX68SJZSsqGhwOHDgGHD3OKC2OVgH7zxPfu\nBZYsAXbvlrsLReevf/kFuHIFEdlmmDp1L5bGn0cv7MFpeGGC0+uYu6y/Kvg/fQqcPy92p9+7JzaN\n3rsHXLmSgXv3qqCgoBosLFLh4mKKdu1s4OQk6pk0agQ4OQGOjiJ7sCTq5uadnWdgyRINKhIWUiqB\ntDTx2wFjzKgYV554YiLw3nuihoc2FM5fAyK6du+O3jt3AkuAVd/nwiL2Npa6vYlRbw3GnTttMWYM\ncOYMEB8vqgS6ugIODkCbNoC9/d+IjT2E/PzJAIDMTFukpwfj7bcDJAVdKdNEZValCgdwxioR/YzE\n8/NF7Q5/f2DmTLlvD6SnA+7uYq7DwQHYv18sQv576927xXT3wYNiWtvL69nD3R2oVq1ocxrtymSM\nsTIynpH47NniRJ4ZM7TT/tKlwGuviYgcFgZYW+POHWDVKvFo2FDMtmzcCNSqVXpzUvLVGWNMF3Qf\nxM+fFyeqnzuneSGm4nK6Z8wATEygVIp1vhUrgBMngGHDxCjc3b18t5GSr64T168DH34odn+amem3\nL4wxvdB9nrinJ3D6tDzFsYvL6TYxQVwc0K0b8OWXwKBBYoEyNLT8ARyQlq+udUlJYjoqMJADOGOV\nmO5H4iYmgL29PG0VTl57ealyugsKRAr3/Pli1mbSpLIP+Ivbkq/VhUgp4uOBnj3FOaPjx+unD4wx\ng2CYpWjL4t49ICBApNSdOAFYW+PKFWDMGMDcXMx9OzuXvTn1aX/BWLJEWgaK1pw9C/TrJ35C8W5M\nxioUKbFT8nRKamoqBg0aBFdXV7Rs2RInT56U2lT5xcYCnTsDo0YBly8jv6Y1vv1WbLIZMeJZ1kl5\nlGlLviGIiAB++okDOGMMgAbTKVOnTkWvXr2wbds25Ofn4+nTp+ovzMsTE9JNmki9VVEHDwJDhuD8\ne+/j8wNP8OSP73H58jto3LgGzp4Vm2+kMJoMlDlz9N0DxpgBkRTE09LScOTIEYSFhYlGqlZF7dq1\n1V88bx7w11/A779L7qQqC8XUFLh0CSc+CcaIX/4pMnJOTw/GpUsBcHKSNvVhsBkojDFWAklB/ObN\nm6hbty5Gjx6Nixcvol27dliyZAksCwtEAQgJCRHz1ps2wScsDD6a9PL53Zf9+yPk0KOXpj4SEko+\nLKI0JRXOYowxbYiOjkZ0dLRGbUha2Dx79ixee+01HD9+HF5eXvjoo49gZWWFr776SjRqYgJ68kTs\nWV+4EBgwoMT2Sj2ooVcvkUbo5QVERaG97yb8738vzwl7e4cgOjqkvF9OkX6Ehu57LgPFT3+LmrGx\nQHCw2KzE2+gZqxQkJYWQBPfv3ycnJyfV348cOUK9e/dW/R0A0dixRGPGlNpWeHgMOTvPIFFDVTyc\nnWdQeHjMs4sUCqLAQCKFgk6cIKpWLb3I9YWPgICZUr4cw5KZSfTFF0R16xKtWEFUUKDvHjHGdERK\nSJaUnVK/fn04OjoiLi4OALB//364ubkVvejkSZGwXYoSs0Ly8sQT1tbAli2IvmCNfv2AL764aXib\nb+Rw4ADQurXIA794UeSAa7qrlTFWoUnOTgkNDcXw4cORm5sLZ2dn/Prrr0UvuHChTPWsi8sKsU5N\nB9zcgGPHgLp1sXu3KHq4eTPwxhut4eWlMJzNN3K4fh0YNw5Ytgzo3VvfvWGMGQm9b/ZRVyHQGgpc\nrNEMDefPASZPxsmTYn/Lzp1Ap05y99aA5OXxFnrGKjGdbvaRy4t1SaohB5HVPZDv9yYweTISE4GB\nA4HVqytIAH/6FPjnH/WvcQBnjJWT3oN4797dsGRJACIc2uJ8bSckVXsFjTwaosn2TcjKAvr3F/VP\n+vTRd081kJMjjksbN04c/bNli757xBirIAziEMbevbsBzlZAzHnxRAM7kEkVjB8PNG0KfP65fvsn\n2cmTwNixQEIC0KIFMHSoOIxZrgJgjLFKzyCCOIBnJ8l7eQGrVmHRIuDKFeDIkdLPrNSZR4/ESRLJ\nyeKRkiL+++qrwIsLu4A4TWjDBhHAzc1131/GWIVnOEF8wwaxvX7FCuw5YY0ffgBOnXoW27UqLU1k\n05w/Lx7Z2SIN5kU5OcC1a4CtLeDiArzying0aKC+XRsb3qjDGNMqvWenvOjaNVGN8I8/gC5dZO7Y\nixQKMfJ/8EDkZ7dtK3aZtmsn/ssYYzokJXbqLIiXurUeQGoq0LEjMG2aWAPUOiLg6lUxqjY1sGqF\njLFKx2CDeHEHLvw2rh5ec7QFhg8HkcgFd3ISx6jJ5v59YO1aUX+ldWsZG2aMMXkZbJ64uq31T+Kn\nomnILNVZmz//LOLt99/LdNP0dGDqVKBlS7EbkhcWGWMVkE6C+Itb602gRBhGIaqeB+Dvj+vXgVmz\ngN9+k2m/y86dYst+ejpw4wbwyy9iyoQxxioYnWSnFB648DOC4II42OEBnsAKy1r44p18YORIEcRd\nXWW4WWqqOIgiLAx44w0ZGmSMMcOl0znxX+KPwQficIeDFs7I2roa5893Q0wMsHevjAX7iAwouZwx\nxspGypy4TkbihVkolu9tBh4Ddy3rIG/1EtSv3w1LlwLnzslccZUDOGOsktBtnnhqqmpDT5a5Ndq2\nBWbPFrvRJbl1C2jYkGtuM8YqBINNMVTno4/EHptNmyTe5OhRUd5w1y6gQweJjTDGmOEw2OmUFx04\nAGzfLg6vkWTHDnHqzfr1HMAZY5WazuchUlOB0aOBVauAOnUkNPDf/wIffABERgL+/rL3jzHGjIn2\ng3h4ODB9uuqvU6eKnZmS4u/q1cCiRaK0Ybt28vWRMcaMlHbnxJOTAXd3UaHQ2xtRUcCECUBsLFCz\npoSGU1LEEWZ2dnJ3mTHG9M7wFjaHDQPq1QMWL8bTp6J0yY8/Aj16yH1HxhgzfoYVxD09xbb3v/4C\nLC0xbRqQlCTWIhljjL3MsLJTLlwA3nwTsLTEuXOikGBsrNbuxhhjlZL2FjabNQO2b0d+vqgNvnCh\nmFkps7//Bj78UGvdY4yxikB7Qfz0acDaGosXixPMRo4sx3uTkoDevXVwtA9jjBk3rS5sJiSIvTin\nTgHOzmV8c0YG4O0NDBgABAfL3TXGGDNYBjUn7u8/E48ff4rp023KHsALCkRGi6cnMGOGtrrGGGMV\nhtaC+L59c1GtWhJatIgF0K3U6wEAP/wAZGaKXZlciZAxxkqltekUQDQbEDALkZFfl+2N2dlAbi5g\nZSV3lxhjzOAZ1HRKoezscpwiX726eDDGGCsTrddOqV69QNu3YIyxSkurQdzZeQYmT/bT5i0YY6xS\n01oQDwiYhSVLeqiOZlPr779FbVrGGGOS6O1kHygUQJs2wE8/AT17yt0FxhgzOoZVAKukZomAwECg\nQQNgyRK5b88YY0bJILNT1PrlF+DGDWDdOr3cnjHGKgrdj8QvXxbb6g8fBlxd5b41Y4wZLSkjcckL\nmwUFBWjTpg369u1bvjfu3g188w0HcMYYk4Hk6ZQlS5agZcuWSE9PL98bP/tM6i0ZY4y9QFIQv3fv\nHnbv3o3g4GB8//33aq8JCQlR/dnHxwc+Pj5SbsUYYxVWdHQ0oqOjNWpD0px4YGAgZsyYgSdPnmDR\nokXYtWtX0UYlzOswxlhlp5M58fDwcNSrVw9t2rThQM0YY3pW7iB+/Phx7Ny5E40bN8bQoUNx8OBB\njCzp2J6EBHEqBGOMMdlplGIYExNT8nSKUgl07w706QN8+qnGnWWMsYpMpymGz9+0WCtWAFlZwEcf\naXobxhhjamhvs8/t20C7dkBMDNCypdy3YIyxCkcvI/FiBQUBH3/MAZwxxrRIe0FcqQSmTdNa84wx\nxrQ5naJU8mHHjDFWDoY1ncIBnDHGtE7rZ2wyxhjTHg7ijDFmxDiIM8aYEeMgzhhjRoyDOGOMGTEO\n4owxZsQ4iDPGmBHjIM4YY0aMg7gR0PT4JvYMf5by4s9T/ziIGwH+H0U+/FnKiz9P/eMgzhhjRoyD\nOGOMGTGtVTFkjDFWfuUNyVUNoROMMcak4ekUxhgzYhzEGWPMiHEQZ4wxIyZ7EI+MjESLFi3QrFkz\nfPvtt3I3X+k4OTnB3d0dbdq0QYcOHfTdHaMzZswY2NnZoXXr1qrnUlJS4OfnBxcXF/j7+yM1NVWP\nPTQe6j7LkJAQODg4oE2bNmjTpg0iIyP12EPjcvfuXbzxxhtwc3NDq1atsHTpUgDl//6UNYgXFBRg\n0qRJiIyMxOXLl7Fx40ZcuXJFzltUOiYmJoiOjsb58+dx+vRpfXfH6IwePfqlwLJgwQL4+fkhLi4O\n3bt3x4IFC/TUO+Oi7rM0MTHBJ598gvPnz+P8+fPo0aOHnnpnfMzMzPDDDz/g0qVLOHnyJJYvX44r\nV66U+/tT1iB++vRpNG3aFE5OTjAzM8OQIUPw559/ynmLSomzfaTr2rUrbGxsijy3c+dOjBo1CgAw\natQo7NixQx9dMzrqPkuAvz+lql+/Pjw9PQEANWvWhKurKxITE8v9/SlrEE9MTISjo6Pq7w4ODkhM\nTJTzFpWOiYkJfH190b59e6xcuVLf3akQHj58CDs7OwCAnZ0dHj58qOceGbfQ0FB4eHhg7NixPDUl\n0a1bt3D+/Hl07Nix3N+fsgZx3uQjv2PHjuH8+fPYs2cPli9fjiNHjui7SxWKiYkJf99qYOLEibh5\n8yYuXLgAe3t7fPrpp/ruktHJyMjAwIEDsWTJEtSqVavIa2X5/pQ1iDdo0AB3795V/f3u3btwcHCQ\n8xaVjr29PQCgbt26ePvtt3leXAZ2dnZ48OABAOD+/fuoV6+enntkvOrVq6cKNOPGjePvz3LKy8vD\nwIEDMWLECPTv3x9A+b8/ZQ3i7du3x/Xr13Hr1i3k5uZi8+bN6Nevn5y3qFQyMzORnp4OAHj69Cmi\noqKKZAYwafr164ewsDAAQFhYmOp/HlZ+9+/fV/35jz/+4O/PciAijB07Fi1btsRHH32ker7c358k\ns927d5OLiws5OzvTN998I3fzlUpCQgJ5eHiQh4cHubm58ecpwZAhQ8je3p7MzMzIwcGBVq9eTcnJ\nydS9e3dq1qwZ+fn5kUKh0Hc3jcKLn+WqVatoxIgR1Lp1a3J3d6e33nqLHjx4oO9uGo0jR46QiYkJ\neXh4kKenJ3l6etKePXvK/f2plQJYjDHGdIN3bDLGmBHjIM4YY0aMgzhjjBkxDuKMMWbEOIgzxpgR\n4yDOGGNG7P8B67fYe1VTegIAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## OLS with dummy variables\n",
      "\n",
      "We generate some artificial data. There are 3 groups which will be modelled using dummy variables. Group 0 is the omitted/benchmark category."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nsample = 50\n",
      "groups = np.zeros(nsample, int)\n",
      "groups[20:40] = 1\n",
      "groups[40:] = 2\n",
      "dummy = (groups[:,None] == np.unique(groups)).astype(float)\n",
      "x = np.linspace(0, 20, nsample)\n",
      "X = np.c_[x, dummy[:,1:], np.ones(nsample)]\n",
      "beta = [1., 3, -3, 10]\n",
      "y_true = np.dot(X, beta)\n",
      "e = np.random.normal(size=nsample)\n",
      "y = y_true + e"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Inspect the data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print X[:5,:]\n",
      "print y[:5]\n",
      "print groups\n",
      "print dummy[:5,:]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.          0.          0.          1.        ]\n",
        " [ 0.40816327  0.          0.          1.        ]\n",
        " [ 0.81632653  0.          0.          1.        ]\n",
        " [ 1.2244898   0.          0.          1.        ]\n",
        " [ 1.63265306  0.          0.          1.        ]]\n",
        "[  9.28223335  10.50481865  11.84389206  10.38508408  12.37941998]\n",
        "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
        " 1 1 1 2 2 2 2 2 2 2 2 2 2]\n",
        "[[ 1.  0.  0.]\n",
        " [ 1.  0.  0.]\n",
        " [ 1.  0.  0.]\n",
        " [ 1.  0.  0.]\n",
        " [ 1.  0.  0.]]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Fit and summary:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "res2 = sm.OLS(y, X).fit()\n",
      "print res.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "                            OLS Regression Results                            \n",
        "==============================================================================\n",
        "Dep. Variable:                      y   R-squared:                       0.933\n",
        "Model:                            OLS   Adj. R-squared:                  0.928\n",
        "Method:                 Least Squares   F-statistic:                     211.8\n",
        "Date:                Sun, 26 Aug 2012   Prob (F-statistic):           6.30e-27\n",
        "Time:                        20:52:53   Log-Likelihood:                -34.438\n",
        "No. Observations:                  50   AIC:                             76.88\n",
        "Df Residuals:                      46   BIC:                             84.52\n",
        "Df Model:                           3                                         \n",
        "==============================================================================\n",
        "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
        "------------------------------------------------------------------------------\n",
        "x1             0.4687      0.026     17.751      0.000         0.416     0.522\n",
        "x2             0.4836      0.104      4.659      0.000         0.275     0.693\n",
        "x3            -0.0174      0.002     -7.507      0.000        -0.022    -0.013\n",
        "const          5.2058      0.171     30.405      0.000         4.861     5.550\n",
        "==============================================================================\n",
        "Omnibus:                        0.655   Durbin-Watson:                   2.896\n",
        "Prob(Omnibus):                  0.721   Jarque-Bera (JB):                0.360\n",
        "Skew:                           0.207   Prob(JB):                        0.835\n",
        "Kurtosis:                       3.026   Cond. No.                         221.\n",
        "==============================================================================\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Draw a plot to compare the true relationship to OLS predictions:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "prstd, iv_l, iv_u = wls_prediction_std(res2)\n",
      "plt.figure()\n",
      "plt.plot(x, y, 'o', x, y_true, 'b-')\n",
      "plt.plot(x, res2.fittedvalues, 'r--.')\n",
      "plt.plot(x, iv_u, 'r--')\n",
      "plt.plot(x, iv_l, 'r--')\n",
      "plt.title('blue: true,   red: OLS')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "<matplotlib.text.Text at 0x46d2bd0>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X18zeX/wPHXYWNmYxNmbseYe6aIymwqGxYqTUVLxKQb\npPQtN1lJJBJLhe4mqagUxkyykR9JptyGuR9bZhtjdnuu3x8fm812trOzc7ad7f18PM7Dds7nc32u\nc868z3Wum/elU0ophBBCWK1q5V0BIYQQpSOBXAghrJwEciGEsHISyIUQwspJIBdCCCsngVwIIayc\nBPIqzs3Nja1btxb6WGRkJM2aNSvjGlVOwcHBBAYGlnc1RCUlgbyK0+l06HS68q4G1apV4+TJk+Vd\nDYsp6WucnJzM+PHjcXV1pXbt2nTp0oWvvvoq3zFFfQi/++67tGrVCkdHR5o1a8YTTzxhatWFFbAp\n7woIkaOotWlZWVnY2FSMP1dL1yUjI4MHH3yQRo0asXv3bpo2bcqvv/7KyJEjSUpK4uWXXwYMfwiH\nhoaycuVKtm7dSsuWLYmPj2f9+vUWq68of9IiF+zZs4eOHTtSr149Ro8eTXp6eqHH3d5qfuaZZ5gx\nY0bu7xs2bMDT0xNnZ2fuu+8+Dhw4YNT1+/TpA0DXrl1xdHRkzZo1REZG0rRpU+bNm4erqyujR48m\nNDQULy8vg3VKT0/n1VdfpUWLFjRq1Ijx48eTlpZWotfCEDc3N+bNm0eXLl1wdHREr9eze/du7r33\nXpydnfH09CQqKir3+FOnTuHt7U2dOnXw9fUlISHB6Gt9/fXXnDt3jjVr1tCiRQuqV6+On58fixcv\n5s033+TatWtFnr937178/Pxo2bIlAC4uLowZM8a0Jy6sggTyKk4pxapVq4iIiCAmJoZjx47xzjvv\nGHVu3hZhdHQ0zz77LMuXLycxMZFx48YxePBgMjMzAXjhhRd44YUXCi1n+/btAPzzzz+kpKQQEBAA\nQHx8PElJSZw9e5Zly5YV2WIHeP311zlx4gR///03J06cIDY2lrffftuo52KM7777jk2bNpGcnMzF\nixd56KGHePPNN0lKSmL+/PkMHTqUy5cvAzB8+HB69OjB5cuXmTFjBqGhoflaz127duW7774r9Dpb\ntmxh4MCB1KpVK9/9jz76KGlpaezatavIevbq1YsVK1Ywf/589u7dS3Z2dimfuajoJJBXcTqdjhdf\nfJEmTZrg7OzMtGnT+Pbbb0tczrJlyxg3bhw9evRAp9Px9NNPU7Nmzdygs2TJEpYsWVKiMqtVq8Zb\nb72Fra0tdnZ2RR6rlGL58uV88MEHODk54eDgwBtvvGEwWJaUTqdjwoQJNGnShJo1a7Jy5UoGDhxI\n//79AXjwwQfp3r07YWFhnD17lr179zJr1ixsbW3x8vJi0KBB+T6I/v77b4P91pcvX8bV1bXA/TY2\nNtSvX7/Y1v2IESMICQlh8+bN+Pj44OLiwrx580rx7EVFVzE6HUW5yjszpXnz5ly4cKHEZZw5c4YV\nK1YQEhKSe19mZiYXL140uV4NGjSgRo0aRh176dIlUlNTueuuu3LvU0qh1+tNvv7t8r5OZ86cYc2a\nNfn6nrOysrj//vu5cOECzs7O+VrULVq04Ny5c0Zdp379+oW+B1lZWSQkJFC/fv1iyxg+fDjDhw8n\nOzubtWvXMmLECDw9PfH19TWqDsK6SItccPbs2Xw/N27cuNDj7O3tSU1Nzf09b5Bu3rw506ZNIykp\nKfd27do1Hn/8cZPrdftAXu3atfNdPy4uLvfn+vXrU6tWLQ4fPpx7/eTkZK5evWry9YuqT/PmzQkM\nDMz3fFNSUnjttddwdXUlKSkpX13PnDlj9MyVBx98kE2bNuU7H+DHH3+kZs2a9OrVy+g6V69encce\ne4wuXbpw6NAho88T1kUCeRWnlGLJkiXExsaSmJjI7NmzDX7l9/T05JtvviE7O5vw8PDcvm2AsWPH\n8umnn7Jnzx6UUly/fp2wsLBiB+ZyuLi4EBMTU+QxXbt25dChQ/z999+kpaURHByc+1i1atUYO3Ys\nkyZN4tKlSwDExsYSERGR75i8dS6Np556ivXr1xMREUF2djZpaWlERkYSGxtLixYt6N69OzNnziQz\nM5Pff/+dDRs2GF12YGAgTZs2JSAggDNnzpCZmcnmzZuZOHEib731Fo6OjrnHZmRkkJaWlnvLysoi\nNDSUjRs3kpKSgl6vZ9OmTRw6dIiePXua5bmLCkiJKs3NzU3NnTtXdejQQTk5OalnnnlG3bhxQyml\n1LZt21SzZs1yj927d6/q2LGjcnR0VIGBgWr48OFqxowZuY+Hh4erHj16KCcnJ+Xq6qqGDRumrl27\nppRS6rnnnlPPPfecwXp8+umnytXVVTk5Oak1a9aoyMjIfNfOMXv2bFW/fn3VvHlztXLlSlWtWjUV\nExOjlFIqLS1NTZ06VbVq1UrVqVNHtW/fXoWEhCillDp79qyqU6eOSkxMNPl12rp1a777/vjjD+Xt\n7a3q1aunGjRooB566CF19uxZpZRSJ0+eVF5eXsrBwUH169dPvfTSSyowMDD33I4dO6pVq1YZvF5i\nYqIaN26ccnFxUbVq1VKdOnVSn3/+eYE66XS6fLcZM2aon376Sd13333K2dlZ1alTR3Xp0kWFhoaa\n9LyFddApJRtLiMrvm2++4fDhw8yePbu8qyKE2RUZyNPS0vD29iY9PZ2MjAyGDBnCnDlzCA4O5rPP\nPqNBgwYAzJkzJ3f0XgghRNkqtkWempqKvb09WVlZ9O7dm/nz57N161YcHR2ZPHlyWdVTCCGEAcUO\ndtrb2wPaoEp2djbOzs5A0cuphRBClKHiOtGzs7NV165dlYODg5oyZYpSSqng4GDVokUL1aVLFzV6\n9GiVlJSU7xxAbnKTm9zkZsLNFMW2yKtVq8b+/fs5f/4827dvJzIykvHjx3Pq1Cn279+Pq6srr7zy\nSoHzlFJyM9Nt5syZ5V6HynST11Nez4p6M5XR88jr1q2Lv78/e/fupWHDhrl5NsaMGcOePXtMroAQ\nQojSKTKQJyQkkJycDMCNGzfYsmUL3bp1y7eibu3atXTu3NmytRRCCGFQkblWLl68yMiRI9Hr9ej1\negIDA3nggQd4+umn2b9/PzqdjpYtW7J06dKyqm+V5OPjU95VqFTk9TQveT3Ln0UWBOl0ulL19wgh\nRFVkauyUXCtCCGHlJJALIYSVk0AuhBBWTgK5EEJYOQnkQghh5SSQCyGElZNALoQQVk4CuRBCWDkJ\n5EIIYeUkkAshhJWTQC6EEFZOArkQQlg5CeRCCGHlJJALIYSVk0AuhBBWTgK5EEJYOQnkQghh5SSQ\nCyGElZNALoQQVk4CuRBCWDkJ5EIIYeUkkAshhJWTQC6EEFZOArkQQlg5CeRCCGHlJJALIYSVk0Au\nhCg7R45AWFh516LSkUAuhLAspWD7dhg8GHx8ICbm1mM3bsCOHeVWtcqiyECelpZGz5498fT0pEOH\nDrzxxhsAJCYm0q9fPzw8PPD19SU5OblMKiuEsCJKwY8/wj33wLPPwsCBcPo0TJhw65gbN8DfXztW\nmEynVNGvYGpqKvb29mRlZdG7d2/mz5/PunXrqF+/Pq+99hrvvfceSUlJzJ0791ahOh3FFCuEqAom\nTYI+fWDIEKhevfBjXFzgr7+gadOyrVsFZGrsLLZrxd7eHoCMjAyys7NxdnZm3bp1jBw5EoCRI0fy\n888/l/jCQojKLSxsO35HHPBZ/A9+A2cSFra98AM7dND6zoXJbIo7QK/Xc+eddxITE8P48ePp2LEj\n8fHxuLi4AODi4kJ8fHyB84KDg3N/9vHxwcfHx2yVFkJUIMeOwb598MQTuXeFhW1n4sTNxMTMzr0v\nJmYaAP7+ffKf3769Fsj79SuT6lYkkZGRREZGlr4gZaTk5GTVs2dP9dtvvyknJ6d8jzk7O+f7vQTF\nCiGskV6v1I4dSg0ZolSDBkrNmZPvYV/faUrr+M5/8/ObXrCsxYuVeu65Mqp4xWZq7Cy2RZ6jbt26\n+Pv789dff+Hi4kJcXByNGjXi4sWLNGzYsPSfKEII6/Dzz/Dee3DpErzyCqxaBTe7YHOkpxceWtLS\nCukn79ULUlMtUdMqo8g+8oSEhNwZKTdu3GDLli1069aNwYMHExoaCkBoaCgPP/yw5WsqhKgYDhyA\nKVPg339h/PgCQRygZs2sQk+1s8sueGePHvC//5m7llVKkbNWDhw4wMiRI9Hr9ej1egIDA5kyZQqJ\niYkMGzaMs2fP4ubmxurVq3FycrpVqMxaEaJKK6yP3N19KosW9S/YRy5ymRo7i51+WJaVEUJUEEeO\nQFQUPPecyUWEhW0nJGQLaWnVsbPL5qWX+kkQL4YEciFE6SilBe/58+HPP7WFO1Ongk5X3jWrMkyN\nnUYPdgohKrGffoLZs+H6dZg8GdasgVq1LHvNoCBt6qK9vTZgmqd7VpSMtMiFEPDJJ9CsmbaMvlrZ\npGDKuNeHGruitF8GDYKnnoJhw8rk2hWVdK0IIazCP9HZLFxcnSdXDsQ3axNnXXrQPHyZtoz/zJny\nrl65stgSfSFEJbFvH7z9dtleMygIfHxQAwaw551wdrgM5eR9T9GmDXT/dxUXegcwuVMEdO4MCQlw\n7VrZ1q+SkEAuRGWm12v5v/v21Vq8tWpp95WR7CP/QlQUuvBwOr4ZQE2/vgw8v5ypU6FeKycyPljC\nyF3jtIRabdpoc9OrguRkbRzCTCSQC1FZrV4NHTvCjBkwZgycPKkt5CmDPvC4OHhzhiL192gAUhu3\nxv6/09y94kVq1HPIPa5ZF2ceTF1Hyn83buVcqcxOndIyQrZqBRs2QHYhC6RMIIFciMqqWjVYskRL\nETtiBNjamv0SYWHb2di0G/ud3NhTvzXL39/G6NFaTE64rOPSpr0QEID9oT/R1b+jwPnVa9oQa+fO\n6S3HtJMOHzZ7HSuEP/+EgABtFWvNmvDPPxAaaji1bwnJ9EMhKqvHHrNo8TmrNz+LrYsn+wE489pi\nsp9248SJltxxB4AH9F9dZDmXG7Yn7fcjMHoAxMZatM7l5vffoXdv+OILcHQ0e/HSIhfCWikF27Zp\nqy/z9HuHhW3Hz286Pj7B+PlNN5wHvJQWL9hC95jOdL0ZxPfQg7F8SXz8FzeDuHEy3NuT8c8RrbVa\nWfM2vfwyYa274ffYexZ5X6RFLoS1yciA77+HDz6AtDR4+WUtkFerVrI84CWRZ/HOf299wr6pa1i+\n7StO0YoX+JiHWUsQy7mCU+EZDotg1609ulU/mV63iiI2Fn74QVsRe9tqWIu9LzlKkTrXIAsVK4QI\nDVWqSROl7r9fqQ0blMrOzvdwifKA32bDhijl6ztNeXvPVL6+09SGDVG3HvT2zi0sEWf1R5vhanTX\nl81yLb+es1WPOw6V9JWoOKKjlXrqKaWcnZWaMEGptLQChxj7vpgaO6VFLoQ1adEC1q+Hbt0KfbhE\necDzMNRi1OvBxqYPdxy0527gQpMe2G1Zz93tXbgUtp2oidMKZDh86aX+Jb6WTpfJjz/+ztChvYs8\nt0LZtg1mzYLjx+GllyAkxGCaAVPfF2NJIBfCmnh7F/lwifKA57F4cQQxMbNZShAeHOMGdkyPeYcn\nn2xO69bw+juruOvXIBp/tiw3WOV0CYSEzMiT4bD4NLU518pLKVsWLPjHugL5f//Bs89qaQWKmRFk\n6vtiNJPa8cWwULFCVH6ZmUp9951SAwYodeNGiU/fsCFKubtPzff13d39jfzdJIXw9p6pQKkd3Jd7\n4imaq65dv1R6valPpuhr3X5r1+5HpQ4eVOrjj817wQqguPclp6vJ1NgpLXIhKoIrV+Czz2DxYnBz\n0zIQmjDv29RWcv2rig+ZyN38AcBBOtCbnfRqtMDsWWwNtU7T0++AGzdg2TJt56GK4O+/4auv4P33\nwcb0cFnU+5K/q2l20QUZYs5PnRwWKlaISmfDhii1sMMAdcXGTm1t1FHtWLjU8hcdO1Ypb2+l7z9A\nRaxOUg8+qNQWW1+10G6Yas9B9R0Bqi5JRrXkTVFY69TFZaXq1euSUikpStWqpVRWltmva7TsbKXC\nwpR64AGlGjdW6t13lUpNtdjlcgZClzJWWuRCWJuclljNmPdZQB3OxzXD/aNpLGqz3aI76WQfPUb1\nHVHoAP2uIAIXr8ZrfTjpW3fQPOQ7PknrQC+7BUa15E1RWOt0VI9ONPlgODhEQP36WhbEVq3Mfu1i\n/fILvP462Nlp34oefxxq1LDoJXMGQj04ZnohZv5wUUpJi1yIQt3W2VyaqYKmiDuapD57OkptsR2g\nFKgrbXsofWKSRa5l0N69Sq1fX+Du9P+S1VUc1I1UvVK+vtrUyvIQGanUtm0F3itLyMpSau1apZyd\nTypQKowBJsdOWdkphKUlJ2vbp3XooPWF32TRKWk308cycCDHVu9nS4eJ1GjfCvd9q2mxcxUEBFBn\ndwQ65zLYlSczU0vgdd998OijEB9/67H0dIiOpkaDutywceTU9nPa61ReOVe8vbXXzVLb2wUFke3l\nw+kOA+nunsScOTB2bBqtWs1gOKtMLla6VoSwlJMnYdEi+PprGDBA+7du3dyHLTklTR07hi5K232n\nxaatJPR5Gb7/B5/OTbUDVhed/ySvsLDtLF4cQXq6DTVrZjFhgq9xXS7Z2dog4ZIlWjfJK6/A4MH5\nBw1TUrQUu0lJxDm3JzHyCO1HjdK+nFiCXg8bN2oDyytXgoND8eeUkKHX6/x5SA8/hvu5KNyAPQ6t\nsDmRgM6mPX369CMkZAGbN5t4UTN/W1BKSdeKEOrDD5W64w6l/vc/pc6dK/QQU6cKFiU1VamlS5WK\nrK11n6Tc0Vylnyz8+sYovI5Tja/j228rtW9f0cfUr69UbKza3f0FFea70OS6FunaNaWWLFHKw0Op\nO+9U6uuvlcrIMLk4Q6tgC3u9mjRZory945SXY7RKqt1Yu7NhQ62b6Tamxk4J5EJYwoUL2gyMYmzY\nEKX8/KZrS9X9phsdIDdsiFJhTTxVdN0W6o873NV3n25WM2Yo1aCBUoMGKbV9XZLSBwQolZRU4DyD\ny/ALUSb9+F5eSv36q/pz1BK1uWWQ+crNsWKF9qE6ZIjWB56n/7ukr0fOOYY+3PLOQNmGtwpjgKpL\nkvrC+RGV3biJUjNnavVIKnxsQgK5EOXh2rUyv2ROINmGd24k+Z5H1cCBsero0eLPK0nr2tDiHW/v\nmdoBV65o3z7mzTP9CY0bp1RIiPpnzw3VqV2m6eUYcuiQUsePF7jb1G8bRX249e49W4HK9958R4B6\nrOckpdLTi62qqbFTBjuFMMW//8Lzz2u5Ty5fLtNLL/owgkYxA2jHUQBiaUwQn5Od/Qlt2xo+r7Cl\n8TExswkJ2WLwHEP9+C31CdpONy1bwv/9H3h5lfyJ5Li5M1CbznYcP2VDRoaJ5RjqV+/QAVq3LnB3\nca+HoXTAhgapDx/245/dYwFIxR7QUvuOYxkpTg4WncYog51CGEsp2LoVFi6EvXth3Dg4eJASJd8u\nqTzpY298vorf39rKu79tpC7fM59XuIddPMsXRqWPNWWWzIQJvsTE3EqMVY1sNtbugM/f8XDfeNi/\nH5o1M/35gZaHPCkJOzto3hxOnNBir9GuXoXPP9dWhO7Yoc1DN0JRr0dRaWdzPtxy8tKkYs8c3mBi\n4hz614mmu9NIhp9cxVKCGMcy6rvPKzaRWGlJIBfCWO+9p810ePllLe90rVoWvVxY2Haa/rCZrkln\nAQhvPobENvcQ3vJePoxZhJ78Abi42S6mzJIpbPGO871jqPnqC2BvX5KnY9i992o3tC1GDx+GDuc2\na8vjX3vN8HknT2opDVasAD8/bSm9kUEcin49DLfWZ/DCC778889KPOKO4YM2M6gvWzjz2GgcPljJ\ngl0HCQlZYPGFVfmY1CFTDAsVK0T5Sk01eaFISQfVNmyIUk2bLlJhaLNP/qCH6thkttqwIcrk2S4l\nPq8MFsXcbupUpd56Sym1ebNSPj6GD1yyRBvAfO01pc6eNelaRb0ehsYGWrcOU+7uSnl4XFGxtvWV\nAhVby0ltXvmzaU/4NqbGziLPOnv2rPLx8VEdOnRQHTt2VIsWLVJKKTVz5kzVpEkT5enpqTw9PdWm\nTZvMUhkhKoQjRww+ZO5ZDkqp3NwnasAAlZ2QqHa9v0PNqqVt2lCXpNzcJ3lni5Rmtkux5+3Zo9SI\nEUqNGWNUmea08mu9evrRFC04u7gYPjA+3iwDzYZeD0OzT1xcDqqdO29+xp07p1QhM4NKwyKB/OLF\niyo6OloppVRKSory8PBQhw8fVsHBwWrBggVmr4wQ5SZnvbS3t0qt30AF+Ewyao5waWc5KKXy7b6T\nVM1Zna7RWr3lMlaB3vBsEXPLzFRq9Wql7r1XqRYtlJo/36wByljHP9qkdjr00yKlg4NSZ86UeR2U\n0t7rJk2W5Jt9sr5WF4skEcvL1NhZZB95o0aNaNSoEQAODg60b9+e2Ju7XCtLrbwSoiylpGg7my9e\nDPXrs8+nP8PP9eLfyLm5h+QMchXVb1pUH2hRg2pxcaA7nIALkFrdgfOvf0TH4MfZ6T8TIgouEzfb\nRgR5ZWdDly7aoO3kyTBkSKlStpZGk75tsbt2mKxsHTbt22ud5jEx0LBhmVw/O1vbgGnhgt54pyTQ\ntdo+0MOxOq7YLnsXP0v3dZvI6Hfr9OnTREdH06tXL3bu3ElISAgrVqyge/fuLFiwAKfbtjgKDg7O\n/dnHxwcfHx9z1VkI8/nsM9i1SxvE7NWLN/rP4N+Tc/MdkhOsTc2NUtgsh+Gs4sSJQbRvDx96vspj\n6d9Qe+MaOt38f3T7bBEwbhs1k1SvDhER0LSp+cs2VkIC/PEHtQYMoJ4uidP/XKX1ihXah0uDBpa9\ndlAQ2UeOce6yPU9e/4yHdb+wPuNDHFs7oBs3HyIi8PjsMzwMbONWGpGRkURGRpa+IGOa7SkpKequ\nu+5Sa9euVUopFR8fr/R6vdLr9WratGlq9OjRZvl6IER5K2oBjKmrHAtbwLOah9XTT59UCQlFn2dK\nP7hBer1Sly+XrgxL+fdfpdzclFJKnajjqbbN+8NsRRc1rnHunFInmnnnvi/Xm3ko/ZAhSkVFlctg\nr6mxs9izMjIylK+vr1q4sPAcCKdOnVKdOnUyS2WEsIi0NG37tNt2nC9MUcHa1NkiqVcy1Be+H6hz\nuCoF6p9a7mrDyjBzPbvipaYqtXy5Uh07lssAplEyM5Wys1Pq+nW1r/1wtW7oV2Yp1tC4xsKFf6nh\nw7WN7w+10GYGqR49yv2DziKBXK/Xq8DAQDVp0qR891+4cCH35w8++EA9+eSTZqmMEGYVF6dUcLA2\n+6FfP6X++6/YU4zZW9FQK/n2/Cc/fvCjCn9gnjpXrZn6p14fdXDqN0r/mHlnORTpwgWlpk/XEjQN\nHKjUli3l0so0WqdOSu3bp/4cNk993c3wZIqSKGz2iROXlU+NTWr+fKWSk5X2fph59kmxEhOV+umn\nAndbJJDv2LFD6XQ61bVr19yphhs3blSBgYGqc+fOqkuXLmrIkCEqLi7OLJURwiwOHlRq5EilnJy0\nPB6HDpXodFO6NArrPknDVm1t/Kg69UPBLHcWl52tVNu2So0fr4pMwFKRBAQo9c03as8epTw9zVNk\nTldZ3vflCg7qiENjo5Kamd3Ro9p74uys1KhRBb4lmho7ixzs7N27N3q9vsD9AwYMKH3nvBCWcvQo\ntGsHCxaYtHze379PiVfiLVq0hZiY2aQyEIDDtOMxfqB55+8IH3pXietQatWqaekDymn2iUnat4fD\nh2k3SEtlk52tjcOa6vx5OHPmfly5QFPOA5CEEyNYif7eXYRbIBe5QVFRWm72PXu01A6HDoGrq9mK\nt6J3WQgjDR1q2fJvy3+ycoMTO3dOAGA4t3JsXMGJhubY7acoKSlw9qw2Te921hTEAXx94fx5HB21\niSqnT4O7ewnOv/m+XMm059XGq/hxqxNeXq3okTyO35L78i8ejGAV9d3nsWiCZXOfFLB3rzatc80a\ny6R2MMe3hdtZqFghbvn3X6VeecWiu5sblGcBT4Kuvhrsn6W6d//MpBktJjt9WqnJk5WqV0/bvKKS\n6d9fqXXrjD8+K0upSx29c1/4I10Ccru8zT77x4JMjZ2SxlZYD6Vgyxbw94fevaFmTW0/yDL078YY\nLu09DcC1mvW4/tUaftlQneDgNri7T8t3rDbvu595K7BrFwwbBnfeqf3+118wd27R51ghY7ftvHYN\nPvoIurZJJfFkEgD67j1oF7WMnGnf/v59CA+fRWRkMOHhsyyXwOrAAZgzxzJlF8PKvnuJKmvTJnj1\nVa3vd9Iky2cfzNN9ol+5is1/OJH4/HQGnPmU4/cEUtfhBA7ff43DzWhRWJZAs2e90+th1iwt09/n\nn4Ojo/nKrmA8WySxPzIZaFnwwaAg0g8c41S8Pc8kLmBKo5VEJy3Dpu/dUKMV1b78EiyweKdQer32\nt7lwofbJ88ILpe/cN4HuZnPevIXqdLKEX5jXn39qzS9L7nCel4+PNkAFbHIMYKr7amYOP87AMY2p\n4Vzb8tevwsLCtrNnwld0OZvAsvu75Nvsee9esPf3ocN/2nujt7Gl2rggmDgR2rQp24p++y289ZaW\nzvfll+Hxx0u9eYSpsVNa5MI69OhRZpeKO5/F1fP2eADH6vbA4etl7HsIdLoyChQxMVqu7X5m7pax\nAjkbOtQ9+QJDGUVExDpOnJjOrl13sH17R86cgW317OE/oGlTqm3fru1SVB5q1oSlS6FPn7JpXBRB\n+shFxZCWpm0McO+9Wt6NshAUpLW8Bw6E5GT+3RjDlo4TSWvhwWe9v+Jq/wA8TkfgNcjJ8v9PldJ2\nt3n0UejZU9tUoQr6ac633BvTlhdYQkcOEcYALp98lY8+qsMLL2ifca12rYKAAK1PuryCOGjvlbd3\nuQdxkK4VUd4uXoRPPtG26fL01L4i+/lpfeGWlqf75FKNxlTPTOeo11jaL3kR505NLH990AL4t99q\nfaxXrmh1RAAaAAAgAElEQVT9/yNHQu2q2X3z3J1jGRu9jxQcc3ffWcNQPurTiaio4LKtjF4PGzbA\njz/Cl1+Wyd+kqbFTWuSi/Cxdqk1PuHwZIiMhPBwGDCiT/zA3bsDZy9pWZem6mpx5ZBK1L53m3qg5\nZRfEQWvN/fknvPmmtpDp+eerbBAHiHd2oi3/koo2kB1PQ15lPrVqWSB9ryHXr8PHH2uLymbNgv79\nDW/sXEFIH7koPw89pA0QWXKGQZ7ZJ6xaRVyaEx9/DIsXZ1DfJphldeJY0s2H0YE96X5HGa70y2vh\nwvK5bgU0ZvIgUqKW8r/seaTgWGDz4rCw7SxeHEF6ug01a2blGwg1i48/1j5U+/TRZgb17l0huk6K\nI10ronLL032yxy0Av+TV3HPPBQ4e/J5z517OPczdfRqLFvlZZo6xUrBzJxw/DqNGmb/8SuZStx4s\nqtaa3x3b3pzG2Q9//z6F7mxv9vdt1y5tE4sSLSk1H5NjpxkWIxVgoWJFJWDKnpemlrnul0j1t727\nUqDSsVUvPrZbJSQYsfWauWRkKPXtt1p61DZtlPriC/OWX1m99JJShWwlWWbvWzkyNXZK14ooM4W1\nqHK2UTO1RZVT5msxl3J335lxNI1Dd/xK3/0bqK30/EU3BrEe++iP6b873eSdfkpk4ULt1qoVTJ+u\ndSOVxQBuZRAYWGifdFHvm9FdLqmpsGIFfPONtkrYzs7ctS8XEshFmTF1z0tjyvTAJ3eWQ6uzMcSd\nc2SWCmYDD6FyxvRvXqtmzcK/upp1P8zq1eGnn6B7d/OVWVUYWDOQs2Xe7VJSYotvIFy8CEuWaLOj\n7r0XZs/W5oFXEtJEEGUmPd0GHXqW8Dw1SM+9vzQt4ZxWWiraDJQ99KA3v/NwnU6sZ/CtIJ7nWhMm\n+Fo+L8qECRLEzczQ+6ZUDQMNhC3aLwsWaLOjkpLg99/h558rxCIec5IWuSgzNWtm0ZW/eYCtZHCr\nNWRKS1ifrdi57CAHDjwFFEwfe4fNjULPs7PLNk9elOxsLSAcPgwzZpS4/qLkDL1v77//W6HH5zYQ\nBg2CZ54xKTe91TBzX71SSgY7ReE2bIhS8+r5qA+ZUKI9L3ONHauyevdRl5p1VQdqdlOna7RRb7y4\nV7VqNaPA1mwzZy4xaX/NYl29qtSiRUq1bKlUr16FbtclylbOIKiObKsfBDU1dsr0Q1GmLnfpxoe2\n7dhx29Sy4ny/NIJ7XnyG5lkXATjX3JOmJ/9CV70aYWHbCQnZkqeVdmu6WmH3m2zWLFi0SJvS+Mor\ncM89ppclzCZi1ToOvbiAwUnn8WQ/13DE3X0qixaZOftkGTA1dkogF2UnJQUaN4a4OKNXLx44AJMn\nx9Fr6yc8p5bThIvsoQfj3Lx456MhZfsf9YcftDzgrVqV3TWrqt9+gyNHtLSwhvz7rzYz6PvvOXv3\nfbxzvSnHbBqZ50O7nEggFxXfunWweDH8+qvhY4KCUMeOcem6PUG1V7HnmBO1a0dw4oQvdUnO1w/u\n5zeD8PBZZVd/UXbCw7VByi1bCn98wQKYNw+ee05La+DiUrb1sxBJYysqPh8fbYNdA25czeRG2C7q\nXThIQ+DDu4NofHo1vr7/x4kTvlzBiSdYnXu8Wed9awVq84v/+EObpibKT/v2WovckBEjtABuyc1F\nrIhMPxRlp06dQpP/x/+bzOYH3yfByR0ungPg3zqu/Dt5JDVqGJ4/bLZ53wkJWv+3m5uW6e7xx81T\nrjBds2aQnAyXLhX+eKNGEsTzkEAuyt7NPODXeviwrd1z1GjfCtsju3mu4WO0Uqf5ngB6Xj3MC9P+\nj7Cw7Zad9/3mm9qHy+nTsHUrbNwIDzxQ+nJF6VSrBm3bQuvWhoO5yCV95KLs5UlklVC/HdW3/coT\nr3xCRMQ7BQ7N6Qc3+wyUHJGRWrrSRo1KX5Ywr337tLnfLVqUd03KjPSRC+thr63CvNSyBw32RYCT\nU7H5T/z9+1hmFoKPj/nLFOZx553lXQOrIV0rwvKysrRkRTfpV67i5xoBZKyPyM1FbrF+8GvXICRE\nS1ol3xJFJSWBXFjejh3a9m03RZ9y4o1Wq2nS8daGEmbvB794EaZN0/Z0jIrSfq5EuTWEyEu6VoTl\nhYfD/ffn/rplS8EN4s2S/yTHW2/Bhx9qU9R27dIGzISoxIoc7Dx37hxPP/00//33HzqdjqCgICZM\nmEBiYiKPP/44Z86cwc3NjdWrV+OUZ7suGewU+Xh6alto3XsvAI96XWLUaw0YNMhC19u3Txsgq8xJ\nkkSlZJGVnXFxccTFxeHp6cm1a9e46667+Pnnn/nyyy+pX78+r732Gu+99x5JSUnMnTu31JURlVBc\nnJZC9L//wMaG1AvJZDVpju7SJRzrV5580EKYg6mxs8g+8kaNGuHp6QmAg4MD7du3JzY2lnXr1jFy\n5EgARo4cyc8//2xClUWVEBGhzcu20Xrxjn2ylcPO95UuiKekaF0n99wDGRlmqqgQ1svoPvLTp08T\nHR1Nz549iY+Px+VmbgMXFxfi4+MLHB8cHJz7s4+PDz4yzatqSk6GRx/N/TXtl3CSe/U3rawLF7Rc\nLZ99pn04LFoENWqYqaJClL3IyEgiIyNLXY5RC4KuXbuGt7c3M2bM4OGHH8bZ2ZmkpKTcx+vVq0di\nYuKtQqVrRRRGKS7WbEHcigi6PdGuZOfOnw/vvqsNYL78smQgFJWSxRYEZWZmMnToUAIDA3n44YcB\nrRUeFxdHo0aNuHjxIg0bNix5jUWVk/D7UbKydHQe2rbkJw8aBKNHQ7165q+YEFauyD5ypRTPPvss\nHTp0YNKkSbn3Dx48mNDQUABCQ0NzA7wQpKUZfCh6czy/t30WG9si5nMbao20bStBXAgDiuxa+f33\n3+nTpw9dunRBd3MxxZw5c7j77rsZNmwYZ8+elemHAvR6CAvTuj/atjWYAnbUKG2D9OefL+TBa9fg\n88+1/u+dO7VMiUJUMbKxhCh76emwcqWW5N/ODqZMgYCA3BkqeSmlZSbdtu22TLZxcdoS+mXLtLwn\nr74KPXuW2VMQoiKRpFmibGVnQ5cu2qBjSIi2crOIJfBHj2rxPd8iy+XL4X//gyefhN27wd3d8vUW\nohKSFrkw3aVL0KCBUYcuXqztv7l8eZ47z5zR9u6sX98y9RPCylhkQZAQAFy/Xvj9RgZxKDy/Ci1a\nSBAXwgwkkIvC5Qxg9u0LN1fxmirzzAXab1ksG+8IYSESyEV+6ena7JFOnWD6dBg7Fr79tlRFnvp4\nEw/Y75IcVkJYiAx2iluUgt69te6Ojz6Cvn0J27iDxQ+9RXq6DTVrZjFhgm+JU8tmrN/M1XsHWqjS\nQggJ5OIWnU7rzL65JiAsbDsTJ24mJmZ27iExMbc2f1i8OKL4AJ+VRfNjv5I080Ptd6XgjTe0nOE1\nJfuhEOYgs1aqqqQkcHYu8hA/v+mFboh8551juXKlYb4A7+4+jUWL/AoE85SIXZwZMI7W1//Bzg5t\n6sojj8CJE2Z5GkJUJjJrRRRPr4dffoE+fWDYsGIPN7Qh8qlTKfmCOEBMzGxCQrYUOPb8F5s53Ly/\nFsRB2y0oz7ZvQojSk66VquDGDQgNhYULtaXvr7wCjz1W7GmGNkSGwrtEcna8ByAoCI4dw+aYntQn\nQm7dv3kzTJhQgsoLIYojgbwqGDRIW3izfDl4eRm9CfGECb7ExEy7rQtlKnXqOJAni3GufDveHzsG\nUVG0Aeofmg2s1uaj//GHNqVRCGE2EsirgvXrCfvtTxbPjiA9/TejZ58Y2hAZYOLEggE+5zEAVcse\nHRBt04Ou395MohUVBd27g6OjeZ+fEFWcBPLKQik4f17LTHWbsN/+NDj7xJhgbugY3XN30jglkQxb\nG5JmfYhfteukrfqJL688yufHVzHHOQiWLqNavZuZMe++W1urL4QwK5m1Yu0yM2HNGi0DYe3aWqv3\ntq4TQ7NP/PxmEB4+i7Cw7cZNJbydj492PUDv4Ei8fSv+dyOY4x17k529llq1YrGzM23uuRBVkWQ/\nrGquXNFydy9apGUNfOstGDiw0P5vQ7NP0tKqFzlXvMjgm5mpfQMArtvWZYz6mvoBD+HV9U/+772F\nJrX+hRCmkUBurUaP1jYeXrsW7rqryEMNzT6xs8tm8eIIA1MJZ9wKvDdnoGBvD6tWoa/jRNgmG5Ky\nnqRpzb388/q3fDzRCWdn8PNbV3x5QgizkkBurVavhurViz8Ow7NPXnqpP++//1uh5+SbSnhzBgrA\niQeC8L+2GgcHHa+8OwuvALjf9tahRbX+hRCWIYG8ItPr4eBBbQOH2xkZxMHw7BN//z4sXhxR6Dl2\ndtna/POVK8m88B+2wD6bHixotIzl/zM8i7Gw1n81sqldI8Po+gohSkYGOyui69dvLeBp3FjbH62a\nZRbh5vSRvxZzCQ+OkYo9C5o155PeabTcvImDtXowJ3k8rzcKxWn1Mlrd6WRUeXlb/080foqP6h3k\njgP7LfIchKgsZLCzMrh4Ucs6uGyZloXwyy/hvvuMXsBjipzWerPAp+mSdAYA31gbNkaNYYzNdvxf\nasvSIHB29i9ReXlb/2/X1XFHa+POF0KUnLTIK5Lnn9da3pMm3ba5pXFMnkYIZPkNxCZiE2dsW/Gi\nxxaenNqKgACwtS3+3GJ17w4ffKDleBFCGGRq7JRAXkkU1qVRICPhbbNPcHLiwgVYsgS++zSZULsg\ndMuWce9AJ/N9CfjvP/Dw0Pb3NMunghCVl2Q/tBY3bmiJo8zM8DTCPBkJc2afbNpEWteejBypbQR0\n9Sps/sOJ3rGruc/fjEEctPzmfftKEBfCgiSQl5X4eAgOBjc3+PhjyDKUWdA0xU77O3AAdTMHeJJt\nAx5NXUmHDhATAyEhJvXkGOfCBRg82EKFCyFABjst7+BBbfbJTz/B449rLeJ27cx+mZxpf0sJyp19\nMpxV2NXMIst/MOk79/JVtWfpUPcA/733FT+PcqJGDbNXo6ApU8rgIkJUbRLILW39emjZEo4f1/bC\ntJCcRT8eMcfwQVu8s9Lei6+dwhi8+Tz2XncxYUrNkmSxFUJYCRnsrETCNkTRYNSz3J0QQ7RtRx6y\n3cajoxswcaIFu07yUkqbA+/gUAYXE6LykcHO8hQfD0uXlt31goK0zIMDB0Jystb/PWo03eYu4e12\ne1lfK4Coqb9z8HwDy/Z/50hLgy++gG7d4M03LXwxIcTtigzko0ePxsXFhc6dO+feFxwcTNOmTenW\nrRvdunUjPDzc4pWssA4c0JJXtWsH0dGQnm7xS4aFbefvHzbnzj654daK615+fLC+NYHXPmH4eCf6\nX1nNpGCn4vZWLr3YWJg+HVq0gB9+gPfeg/nzLXxRIcTtigzko0aNKhCodTodkydPJjo6mujoaPr3\n72/g7EosMhL69dM2EW7dWtsR/tNPoWbhe1kaEha2HT+/6fj4BOPnN52wsO3FHj9x4mZikzoAcJ1a\nTL4ykwFt/+TutVP5NfoOhg8vo5l+N25Az57aN4Lt22HjRu31sFAqASGEYUUOdnp5eXH69OkC91f5\n/u/YWHj6aW0WiolTP4rNA17I4p2cueLDSWYlT/EUX3MFZ/ycZ+DlNcssT81otWrByZMmP38hhPmY\nNGslJCSEFStW0L17dxYsWICTU8FESsHBwbk/+/j44OPjY2odK54RI0pdRLF5wPOkjlVjg9jw9Gr2\n7HkWgCs4MYgNuedZNEXspUuQmAht2xZ8TIK4EKUSGRlJZGRk6QtSxTh16pTq1KlT7u/x8fFKr9cr\nvV6vpk2bpkaPHl3gHCOKrfj27lVq0iSlsrIsUry390ylTfPIf/P2nqlUdrZS3bsrBSrNrq66yz1J\n3XmnUp07f1/oOX5+081fwb//Vmr0aKWcnJT64APzly+EKMDU2FniDs2GDRui0+nQ6XSMGTOGPXv2\nlP7TpKLIztZ23OnTBx5+GFxdzb4CM0feBTzb8CGMgbhwkUcu/Ulmu05cPJPOXzZ3M9bnGB984cTe\nvTBnTiPc3aflK0fbIKKfeSqVnQ3r1sH998OAAdCqlfbN4OWXzVO+EMIiSty1cvHiRVxdXQFYu3Zt\nvhktVu3777k+aTJnb8DqxvfwR7v7eKFjL/xLOIBprMIW8BzWtebf+Pt4NG0Jbs/4MHGSjhV5pg4W\ntUGEWWRmaoO2Y8bAY49J14kQVqLIQP7kk08SFRVFQkICzZo146233iIyMpL9+/ej0+lo2bIlS8ty\n/rQF7ToVx4Lq/fjxyldwBTgCx06VftNgQ6llc8q0H/k9XIZ91Tszss7PPDWlFSuCMDh1MO+5Zmdn\np80+EUJYFVnZeZOf33QiIt4p5P4ZhIebNiOksN13prdox/yZveg5bBgrVsDnC5J593IQV+Yt4+Fn\nnCw/dVApbfpkVpY2hVIIUWHIys7iZGTAypXg6wupqQUetsSmwTkzUzzQuk8GsonIM5+hf3ERbm4Q\nEQELv3SiX+JqAsZaOIinpWk7Dnl6ahtYXL1qwYsJIcpS5U+alZCgLZ//+GNo317bfcfOrsBhhW0a\nDDc3ITZReroNLsTRhPMAJFGXQL4mum4Ddm0vo/wnaWkwZ47W933nnTBvntYSl4U7QlQalTuQf/QR\nzJgBjz4KmzYVvhv9TTmDj/l32JnKSy+ZvnK1Ro1svuIZttOH47RhON9yBSf8usygdeteBs8rzZZt\nhVQC9HqtO6V9e9PKEEJUaJU7kPv6wrBh0LBhsYeWekZInpWYqZ+tIvQXJw4ffoMhtolkZDbNPay4\nD4diV3yWVLVqMKuMV30KIcpU5RjszMoCm3L8TLpxQ9v1ft8+AH6uEcBXA1YzeTJcvbqdjz7akufD\noV+RAdmkQderV7X+71q1tA8UIYRVMjV2WneL/MwZrfvk22/h6NGyz4MdFwcff0zmR5+SkmZDPeBM\nwx502bSMn+/MOagPDz1kfEu6RIOuJ09q+7StWKH1e7/6aomfghDC+llfIFcKdu6EDz+EbdvgmWdg\nxw6jg3ip+p9zuk+yslDNW5C5biMb6z7Jwho7eGSCC8/vD6LFV8ugkNwzxjJq0DUtDYYP17IOjhkD\n+/dDs2YmX1MIYeXMlCIgHwsVq5kxQ6k2bZQKCVHq6tUSnbphQ5Ryd5+aL0+Ju/tUtWFDlHEFeHvn\nnnikRid1f9cE9c03SmVklPxplKyObxSs45o1Sl27Zr4LCyHKnamx0/r6yK9cAUdHk6bPlWbRz4UL\nkOI1kLYnN3HcqQeXvongngFOFtn/MixsOyEhxverCyEqh8rXR37ypJa06XZ165pcpFH9zzndJ9Wr\nw8iRRHd+moULYcMGGBuwiunngmizahltStF9Uhx//z74t6yvdR/VqwcSxIUQRahYq0Jysg96e0Pf\nvnDtmlmLN6r/OTpaywP+22+ce24WgwcpOnWCmBh4b6kTjhtXl6oPvEhKacs9+/fXMhA2aQKTJ1vm\nWkKISqNidK0kJ8Pnn2szUFxdYeJEbRFPMWvWixq4LOwxoMAcbXf3qSxa1B9/2zSy576P2rEdm6wM\njtbqxoFFv5VN/hPQUgjcfbe2eOfll+HJJwtdgSqEqLxM7VqpGIH8zTe1rpSJE6FHD6NOKWzhjLv7\nNBYt8gMKC9i3HtM9N4nGKYlk2NqQ9NEiOnv5c2DEXH7Z24TrXn68f+NFGvy0DJ2zcS1vs63E/Ocf\n6NwZi3S8CyEqPJPHF80w0FqAhYrNx9d3msHdcop6TCmVb/bJ7uYBytlZqRdfVOr48ZLXw6SZMOnp\nJj1nIUTlZmrsLLs+8pQUbeGOmb4AFDVwaegx14QE9PMXEH/NHoD9tj3YPXoZMTHaupqiklgZ2vHe\n8N6bW/IXoNdrub4ffBBefNHIZymEEMWz/KyVkye1vu/QUG0Ac/BgqF073yGmdE0UNXCpbn5YLCUI\nD45hz3WuU5uOf//Dh6df5ZfG37Ci5zg6rl+GZ4Piu0+Kyn9S7EyY1FT4+mttBoqdnTZ4+fjjxV5T\nCCGMZt4vBhpAqZ07lRo8WKk77lBqyhSlTp8u9FhTF+kUtXAm57FDtL/VhUJ31fvuWBUVpZReX7Ln\nY3I3TmamUi1bKjVokFK//VbyCwshqhRTQ7LlWuQxMTBwIKxaVaAFnpfhrokZRbbKi8tWePJkbfST\nf4As2G/bnqPvz2XHxMYmPZWiWt1TptxvOP2tjQ3s3avNBRdCCAuxXCAPDDTqsNLszOPv3wf/X1Zq\nC3iq2aO/5xXWr4eFC+HYsbvInvYHHvuD8PxqGZ4Wyn/i798H9HpWLJxMvL5OwfS3EsSFEBZW7is7\nS7Uzj14Pu3fDgQMAhLcIIthjNa+8AgEBYGvrBKwudR0NbToxMagvLFuG/8KF+Ht5wbIPSn0tIYQo\nqXIP5CbtzHPjBnz9NZnvLyT7bCx2wHGn7tRduYy9A80/Dfv2bhxX3RXebpRIm+dGQK9e8Mkn2mpU\nIYQoBxViQVCxSaLy7L7D/feT+e48DtTuRXDyy7R7wpMZseNwXFW69LFG0+u1RTve3toCprZtLX9N\nIUSVYN0rO4vj46PlPwGinXyYVPNT/Ce3ZexYcHY232WMVt47EgkhKqXKl/3wpuvXIeGSPS2AQ/Y9\nODFvLb+WRf6T9HRtByIPj4KPSRAXQlQgFSf7YVCQ1vL284MPPyStTz9mvJ6JmxtMbbmKS30D6HA+\ngoCxFg7ily/DO++AmxssWmTBCwkhhHmUWdOy2NWbBw/Crl0AJO44wNPVV9Gyiw27dkHr1uaZfVKk\nY8e0eYvffadlXtyyBTp1suw1hRDCDMokkBe1xN3fvw/6ue/BH39SDThi05nNr23n64lOODtr577w\nghkyCxZFKRg3Dnr3hiNHoFEj85YvhBAWVCaBPGf1Zk7uk1TsGR6zioULF3P2bB/WfzqI6u0D+Mjh\ndVqvX0b7m/lPivsAMBudDn77TdLHCiGsUpkE8pzVmx4cwwdt9slSggjc/iWOjvD6ig54eYFOl7/7\nxNTl+wZdvaq1uHv2LPiYBHEhhJUqcrBz9OjRuLi40Llz59z7EhMT6devHx4eHvj6+pKcnFzsRRqR\nzNvMoDt/ArCHHoxjGb16LWPtWujTp/A4Wprl+/mcOwdTpkDLlrByZcnOFUKICq7IQD5q1CjCw8Pz\n3Td37lz69evHsWPHeOCBB5g7d27hJw8cCL//jnpmFF/+8SUNq5/Cmyi+JwBfIqjvPo///e+uIitX\nquX7APv2wYgR0LWrth/oX39piceFEKISKXZB0OnTpxk0aBAHbuYzadeuHVFRUbi4uBAXF4ePjw9H\njx7NX6hOhwIybexY7PQmYY2DuMv3P/7++zsyMnSFr94sROHbud3cY9OYrpWnntKC+NixZbPqUwgh\nSqHMFgTFx8fj4uICgIuLC/Hx8YUe93y1Rux0e4o7e6cz4+kD9O3rA7xVomsVl6q2WNKNIoSowCIj\nI4mMjCx1OSVukTs7O5OUlJT7eL169UhMTMxfqE5H3NEkXNqWQSv48mX480/oX0SSLSGEsAKmtshL\nvLIzp0sF4OLFizRs2LDQ456eMD93X0uLiInR9r5s0wbCwix3HSGEqOBKHMgHDx5MaGgoAKGhoTz8\n8MOFHhcR8Q4TJ242fzDftQuGDtXSx9atC4cOyQCmEKJKK7Jr5cknnyQqKoqEhARcXFx4++23GTJk\nCMOGDePs2bO4ubmxevVqnG4bSNTpdIBWrJ/fDMLDZ5mvxlOmQIsWMGpUkVvICSGEtalwaWxzArm3\ndzCRkcHmvoQQQlQ6ZdZHXlJGz/nO69IlWG3hJFlCCFFJWDSQa1u29TP+hOPHYfx4bdedyEgtmZUQ\nQogiWSyQ+/nNMH7hzu7d8MgjcN990KABHD0KH38s+U+EEMIIFWOrt/feAwcHeOYZGcAUQlRZFW6w\n0wLFCiFEpVZhBztzJSTAF1+U2eWEEKKqsHwgP3lSW4Hp4QH/93+QkWHxSwohRFViuUC+bx8MG6Zt\n4lC3Lhw+DJ99BjVqWOySQghRFVluh6C//tJmoXzxhTaQKYQQwiJksFMIISqIij/YKYQQwiIkkAsh\nhJWTQC6EEFZOArkQQlg5CeRCCGHlJJALIYSVk0AuhBBWTgK5EEJYOQnkQghh5SSQCyGElZNALoQQ\nVk4CuRBCWDkJ5EIIYeUkkAshhJWTQC6EEFZOArkQQlg5CeRCCGHlJJBbgcjIyPKuQqUir6d5yetZ\n/kwO5G5ubnTp0oVu3bpx9913m7NO4jbyH8W85PU0L3k9y5/Jmy/rdDoiIyOpV6+eOesjhBCihErV\ntSIbLAshRPnTKROjcatWrahbty7Vq1dn3LhxjB079lahOp3ZKiiEEFWJKSHZ5K6VnTt34urqyqVL\nl+jXrx/t2rXDy8vL5IoIIYQwjcldK66urgA0aNCARx55hD179pitUkIIIYxnUiBPTU0lJSUFgOvX\nrxMREUHnzp3NWjEhhBDGMalrJT4+nkceeQSArKwsRowYga+vr1krJoQQwjgmtchbtmzJ/v372b9/\nPwcPHuSNN97IfSw8PJx27drRpk0b3nvvPbNVtKqS+fqlM3r0aFxcXPJ9Y0xMTKRfv354eHjg6+tL\ncnJyOdbQuhT2egYHB9O0aVO6detGt27dCA8PL8caWo9z587Rt29fOnbsSKdOnVi8eDFg2t+nWVd2\nZmdn8+KLLxIeHs7hw4f59ttvOXLkiDkvUeXkzNePjo6WcQgTjBo1qkBgmTt3Lv369ePYsWM88MAD\nzJ07t5xqZ30Kez11Oh2TJ08mOjqa6Oho+vfvX061sy62trYsXLiQQ4cOsXv3bpYsWcKRI0dM+vs0\nayDfs2cPrVu3xs3NDVtbW5544gl++eUXc16iSpJZQKbz8vLC2dk5333r1q1j5MiRAIwcOZKff/65\nPKpmlQp7PUH+Rk3RqFEjPD09AXBwcKB9+/bExsaa9Pdp1kAeGxtLs2bNcn9v2rQpsbGx5rxElaPT\n6SbOzs4AAAHrSURBVHjwwQfp3r07y5cvL+/qVArx8fG4uLgA4OLiQnx8fDnXyPqFhITQtWtXnn32\nWemqMsHp06eJjo6mZ8+eJv19mjWQy0Ig89u5cyfR0dFs2rSJJUuWsGPHjvKuUqWi0+nk77aUxo8f\nz6lTp9i/fz+urq688sor5V0lq3Lt2jWGDh3KokWLcHR0zPeYsX+fZg3kTZo04dy5c7m/nzt3jqZN\nm5rzElWOzNc3PxcXF+Li4gC4ePEiDRs2LOcaWbeGDRvmBpwxY8bI32gJZGZmMnToUAIDA3n44YcB\n0/4+zRrIu3fvzvHjxzl9+jQZGRl8//33DB482JyXqFJkvr5lDB48mNDQUABCQ0Nz/wMJ01y8eDH3\n57Vr18rfqJGUUjz77LN06NCBSZMm5d5v0t+nMrONGzcqDw8P5e7urt59911zF1+lnDx5UnXt2lV1\n7dpVdezYUV5PEzzxxBPK1dVV2draqqZNm6ovvvhCXb58WT3wwAOqTZs2ql+/fiopKam8q2k1bn89\nP//8cxUYGKg6d+6sunTpooYMGaLi4uLKu5pWYceOHUqn06muXbsqT09P5enpqTZt2mTS36fJSbOE\nEEJUDLJDkBBCWDkJ5EIIYeUkkAshhJWTQC6EEFZOArkQQlg5CeRCCGHl/h96a7tFOscRjAAAAABJ\nRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Joint hypothesis test\n",
      "\n",
      "### F test\n",
      "\n",
      "We want to test the hypothesis that both coefficients on the dummy variables are equal to zero, that is, $R \\times \\beta = 0$. An F test leads us to strongly reject the null hypothesis of identical constant in the 3 groups:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "R = [[0, 1, 0, 0], [0, 0, 1, 0]]\n",
      "print np.array(R)\n",
      "print res2.f_test(R)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[0 1 0 0]\n",
        " [0 0 1 0]]\n",
        "<F test: F=array([[ 145.49268198]]), p=[[ 0.]], df_denom=46, df_num=2>\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Small group effects\n",
      "\n",
      "If we generate artificial data with smaller group effects, the T test can no longer reject the Null hypothesis: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "beta = [1., 0.3, -0.0, 10]\n",
      "y_true = np.dot(X, beta)\n",
      "y = y_true + np.random.normal(size=nsample)\n",
      "res3 = sm.OLS(y, X).fit()\n",
      "print res3.f_test(R)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<F test: F=array([[ 1.22491119]]), p=[[ 0.30318644]], df_denom=46, df_num=2>\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Multicollinearity\n",
      "\n",
      "The Longley dataset is well known to have high multicollinearity, that is, the exogenous predictors are highly correlated. This is problematic because it can affect the stability of our coefficient estimates as we make minor changes to model specification. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from statsmodels.datasets.longley import load\n",
      "y = load().endog\n",
      "X = load().exog\n",
      "X = sm.tools.add_constant(X)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Fit and summary:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ols_model = sm.OLS(y, X)\n",
      "ols_results = ols_model.fit()\n",
      "print ols_results.summary()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "                            OLS Regression Results                            \n",
        "==============================================================================\n",
        "Dep. Variable:                      y   R-squared:                       0.995\n",
        "Model:                            OLS   Adj. R-squared:                  0.992\n",
        "Method:                 Least Squares   F-statistic:                     330.3\n",
        "Date:                Sun, 26 Aug 2012   Prob (F-statistic):           4.98e-10\n",
        "Time:                        20:52:53   Log-Likelihood:                -109.62\n",
        "No. Observations:                  16   AIC:                             233.2\n",
        "Df Residuals:                       9   BIC:                             238.6\n",
        "Df Model:                           6                                         \n",
        "==============================================================================\n",
        "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
        "------------------------------------------------------------------------------\n",
        "x1            15.0619     84.915      0.177      0.863      -177.029   207.153\n",
        "x2            -0.0358      0.033     -1.070      0.313        -0.112     0.040\n",
        "x3            -2.0202      0.488     -4.136      0.003        -3.125    -0.915\n",
        "x4            -1.0332      0.214     -4.822      0.001        -1.518    -0.549\n",
        "x5            -0.0511      0.226     -0.226      0.826        -0.563     0.460\n",
        "x6          1829.1515    455.478      4.016      0.003       798.788  2859.515\n",
        "const      -3.482e+06    8.9e+05     -3.911      0.004      -5.5e+06 -1.47e+06\n",
        "==============================================================================\n",
        "Omnibus:                        0.749   Durbin-Watson:                   2.559\n",
        "Prob(Omnibus):                  0.688   Jarque-Bera (JB):                0.684\n",
        "Skew:                           0.420   Prob(JB):                        0.710\n",
        "Kurtosis:                       2.434   Cond. No.                     4.86e+09\n",
        "==============================================================================\n",
        "\n",
        "The condition number is large, 4.86e+09. This might indicate that there are\n",
        "strong multicollinearity or other numerical problems.\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/local/lib/python2.7/dist-packages/scipy/stats/stats.py:1199: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
        "  int(n))\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Condition number\n",
      "\n",
      "One way to assess multicollinearity is to compute the condition number. Values over 20 are worrisome (see Greene 4.9). The first step is to normalize the independent variables to have unit length: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "norm_x = np.ones_like(X)\n",
      "for i in range(int(ols_model.df_model)):\n",
      "    norm_x[:,i] = X[:,i]/np.linalg.norm(X[:,i])\n",
      "norm_xtx = np.dot(norm_x.T,norm_x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Then, we take the square root of the ratio of the biggest to the smallest eigen values. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "eigs = np.linalg.eigvals(norm_xtx)\n",
      "condition_number = np.sqrt(eigs.max() / eigs.min())\n",
      "print condition_number"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "56240.8699497\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Dropping an observation\n",
      "\n",
      "Greene also points out that dropping a single observation can have a dramatic effect on the coefficient estimates: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ols_results2 = sm.OLS(y[:-1], X[:-1,:]).fit()\n",
      "print \"Percentage change %4.2f%%\\n\"*7 % tuple([i for i in ols_results.params/ols_results2.params*100 - 100])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Percentage change -173.43%\n",
        "Percentage change 31.04%\n",
        "Percentage change 3.48%\n",
        "Percentage change 7.83%\n",
        "Percentage change -199.54%\n",
        "Percentage change 15.39%\n",
        "Percentage change 15.40%\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 21
    }
   ],
   "metadata": {}
  }
 ]
}