{
 "metadata": {
  "name": "",
  "signature": "sha256:5db438f961adadf8fe26821d384aca07c66313bbe5deae4565558e80a8ab58ad"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Ordinary Least Squares"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import print_function\n",
      "import numpy as np\n",
      "import statsmodels.api as sm\n",
      "import matplotlib.pyplot as plt\n",
      "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
      "\n",
      "np.random.seed(9876789)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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": []
    },
    {
     "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": []
    },
    {
     "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": []
    },
    {
     "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": []
    },
    {
     "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.column_stack((x, np.sin(x), (x-5)**2, np.ones(nsample)))\n",
      "beta = [0.5, 0.5, -0.02, 5.]\n",
      "\n",
      "y_true = np.dot(X, beta)\n",
      "y = y_true + sig * np.random.normal(size=nsample)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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": []
    },
    {
     "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": []
    },
    {
     "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": [
      "prstd, iv_l, iv_u = wls_prediction_std(res)\n",
      "\n",
      "fig, ax = plt.subplots(figsize=(8,6))\n",
      "\n",
      "ax.plot(x, y, 'o', label=\"data\")\n",
      "ax.plot(x, y_true, 'b-', label=\"True\")\n",
      "ax.plot(x, res.fittedvalues, 'r--.', label=\"OLS\")\n",
      "ax.plot(x, iv_u, 'r--')\n",
      "ax.plot(x, iv_l, 'r--')\n",
      "ax.legend(loc='best');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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",
      "\n",
      "dummy = sm.categorical(groups, drop=True)\n",
      "x = np.linspace(0, 20, nsample)\n",
      "# drop reference category\n",
      "X = np.column_stack((x, dummy[:,1:]))\n",
      "X = sm.add_constant(X, prepend=False)\n",
      "\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": []
    },
    {
     "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": []
    },
    {
     "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": []
    },
    {
     "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",
      "\n",
      "fig, ax = plt.subplots(figsize=(8,6))\n",
      "\n",
      "ax.plot(x, y, 'o', label=\"Data\")\n",
      "ax.plot(x, y_true, 'b-', label=\"True\")\n",
      "ax.plot(x, res2.fittedvalues, 'r--.', label=\"Predicted\")\n",
      "ax.plot(x, iv_u, 'r--')\n",
      "ax.plot(x, iv_l, 'r--')\n",
      "legend = ax.legend(loc=\"best\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can also use formula-like syntax to test hypotheses"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(res2.f_test(\"x2 = x3 = 0\"))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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",
      "\n",
      "res3 = sm.OLS(y, X).fit()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(res3.f_test(R))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(res3.f_test(\"x2 = x3 = 0\"))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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_pandas\n",
      "y = load_pandas().endog\n",
      "X = load_pandas().exog\n",
      "X = sm.add_constant(X)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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": []
    },
    {
     "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 = X.values\n",
      "for i, name in enumerate(X):\n",
      "    if name == \"const\":\n",
      "        continue\n",
      "    norm_x[:,i] = X[name]/np.linalg.norm(X[name])\n",
      "norm_xtx = np.dot(norm_x.T,norm_x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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": []
    },
    {
     "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.ix[:14], X.ix[:14]).fit()\n",
      "print(\"Percentage change %4.2f%%\\n\"*7 % tuple([i for i in (ols_results2.params - ols_results.params)/ols_results.params*100]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can also look at formal statistics for this such as the DFBETAS -- a standardized measure of how much each coefficient changes when that observation is left out."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "infl = ols_results.get_influence()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In general we may consider DBETAS in absolute value greater than $2/\\sqrt{N}$ to be influential observations"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "2./len(X)**.5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(infl.summary_frame().filter(regex=\"dfb\"))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}