{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantile regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This example page shows how to use ``statsmodels``' ``QuantReg`` class to replicate parts of the analysis published in \n",
    "\n",
    "* Koenker, Roger and Kevin F. Hallock. \"Quantile Regressioin\". Journal of Economic Perspectives, Volume 15, Number 4, Fall 2001, Pages 143–156\n",
    "\n",
    "We are interested in the relationship between income and expenditures on food for a sample of working class Belgian households in 1857 (the Engel data). \n",
    "\n",
    "## Setup\n",
    "\n",
    "We first need to load some modules and to retrieve the data. Conveniently, the Engel dataset is shipped with ``statsmodels``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import print_function\n",
    "import patsy\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.regression.quantile_regression import QuantReg\n",
    "\n",
    "data = sm.datasets.engel.load_pandas().data\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Least Absolute Deviation\n",
    "\n",
    "The LAD model is a special case of quantile regression where q=0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mod = smf.quantreg('foodexp ~ income', data)\n",
    "res = mod.fit(q=.5)\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the results\n",
    "\n",
    "We estimate the quantile regression model for many quantiles between .05 and .95, and compare best fit line from each of these models to Ordinary Least Squares results. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare data for plotting\n",
    "\n",
    "For convenience, we place the quantile regression results in a Pandas DataFrame, and the OLS results in a dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "quantiles = np.arange(.05, .96, .1)\n",
    "def fit_model(q):\n",
    "    res = mod.fit(q=q)\n",
    "    return [q, res.params['Intercept'], res.params['income']] + \\\n",
    "            res.conf_int().loc['income'].tolist()\n",
    "    \n",
    "models = [fit_model(x) for x in quantiles]\n",
    "models = pd.DataFrame(models, columns=['q', 'a', 'b','lb','ub'])\n",
    "\n",
    "ols = smf.ols('foodexp ~ income', data).fit()\n",
    "ols_ci = ols.conf_int().loc['income'].tolist()\n",
    "ols = dict(a = ols.params['Intercept'],\n",
    "           b = ols.params['income'],\n",
    "           lb = ols_ci[0],\n",
    "           ub = ols_ci[1])\n",
    "\n",
    "print(models)\n",
    "print(ols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First plot\n",
    "\n",
    "This plot compares best fit lines for 10 quantile regression models to the least squares fit. As Koenker and Hallock (2001) point out, we see that:\n",
    "\n",
    "1. Food expenditure increases with income\n",
    "2. The *dispersion* of food expenditure increases with income\n",
    "3. The least squares estimates fit low income observations quite poorly (i.e. the OLS line passes over most low income households)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.arange(data.income.min(), data.income.max(), 50)\n",
    "get_y = lambda a, b: a + b * x\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "for i in range(models.shape[0]):\n",
    "    y = get_y(models.a[i], models.b[i])\n",
    "    ax.plot(x, y, linestyle='dotted', color='grey')\n",
    "    \n",
    "y = get_y(ols['a'], ols['b'])\n",
    "\n",
    "ax.plot(x, y, color='red', label='OLS')\n",
    "ax.scatter(data.income, data.foodexp, alpha=.2)\n",
    "ax.set_xlim((240, 3000))\n",
    "ax.set_ylim((240, 2000))\n",
    "legend = ax.legend()\n",
    "ax.set_xlabel('Income', fontsize=16)\n",
    "ax.set_ylabel('Food expenditure', fontsize=16);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second plot\n",
    "\n",
    "The dotted black lines form 95% point-wise confidence band around 10 quantile regression estimates (solid black line). The red lines represent OLS regression results along with their 95% confindence interval.\n",
    "\n",
    "In most cases, the quantile regression point estimates lie outside the OLS confidence interval, which suggests that the effect of income on food expenditure may not be constant across the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = models.shape[0]\n",
    "p1 = plt.plot(models.q, models.b, color='black', label='Quantile Reg.')\n",
    "p2 = plt.plot(models.q, models.ub, linestyle='dotted', color='black')\n",
    "p3 = plt.plot(models.q, models.lb, linestyle='dotted', color='black')\n",
    "p4 = plt.plot(models.q, [ols['b']] * n, color='red', label='OLS')\n",
    "p5 = plt.plot(models.q, [ols['lb']] * n, linestyle='dotted', color='red')\n",
    "p6 = plt.plot(models.q, [ols['ub']] * n, linestyle='dotted', color='red')\n",
    "plt.ylabel(r'$\\beta_{income}$')\n",
    "plt.xlabel('Quantiles of the conditional food expenditure distribution')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}