<jupyter_start><jupyter_text>Quantile regression SetupWe first load some modules and retrieve the data. Conveniently, the Engel dataset is shipped with statsmodels.<jupyter_code>%matplotlib inline

from __future__ import print_function
import patsy
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
from statsmodels.regression.quantile_regression import QuantReg
data = sm.datasets.engel.load_pandas().data
data.head(10)<jupyter_output><empty_output><jupyter_text>Least Absolute DeviationThe LAD model is a special case of quantile regression where q=0.5<jupyter_code>mod = smf.quantreg('foodexp ~ income', data)
res = mod.fit(q=.5)
print(res.summary())<jupyter_output>QuantReg Regression Results                          
==============================================================================
Dep. Variable:                foodexp   Pseudo R-squared:               0.6206
Model:                       QuantReg   Bandwidth:                       64.51
Method:                 Least Squares   Sparsity:                        209.3
Date:                Tue, 28 Aug 2018   No. Observations:                  235
Time:                        15:58:41   Df Residuals:                      233
                                        Df Model:                            1
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     81.4823     14.634      5.568      0.000      52.649     110.315
income         0.5602      0.013     42.516      0.0[...]<jupyter_text>Visualizing the resultsWe 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.  Prepare data for plottingFor convenience, we place the quantile regression results in a Pandas DataFrame, and the OLS results in a dictionary.<jupyter_code>quantiles = np.arange(.05, .96, .1)

def fit_model(q):
    res = mod.fit(q=q)
    return [q, res.params['Intercept'], res.params['income']] + res.conf_int().loc['income'].tolist()

models = [fit_model(x) for x in quantiles]
models = pd.DataFrame(models, columns = ['q', 'a', 'b', 'lb', 'ub'])
ols = smf.ols('foodexp ~ income', data).fit()
ols_ci = ols.conf_int().loc['income'].tolist()
ols = dict(a = ols.params['Intercept'], 
          b = ols.params['income'],
          lb = ols_ci[0], 
          ub = ols_ci[1])
print(models)
print(ols)<jupyter_output>q           a         b        lb        ub
0  0.05  124.880098  0.343361  0.268632  0.418090
1  0.15  111.693660  0.423708  0.382780  0.464636
2  0.25   95.483539  0.474103  0.439900  0.508306
3  0.35  105.841294  0.488901  0.457759  0.520043
4  0.45   81.083647  0.552428  0.525021  0.579835
5  0.55   89.661370  0.565601  0.540955  0.590247
6  0.65   74.033435  0.604576  0.582169  0.626982
7  0.75   62.396584  0.644014  0.622411  0.665617
8  0.85   52.272216  0.677603  0.657383  0.697823
9  0.95   64.103964  0.709069  0.687831  0.730306
{'a': 147.47538852370573, 'b': 0.48517842367692354, 'lb': 0.4568738130184233, 'ub': 0.5134830343354237}<jupyter_text>First plotThis 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:1). Food expenditure increases with income.2). The dispersion of food expenditure increases with income.3). The least squares estimates fit low income observations quite poorly (i.e. the OLS line passes over most low income households).<jupyter_code>x = np.arange(data.income.min(), data.income.max(), 50)
get_y = lambda a, b: a + b * x

fig, ax = plt.subplots(figsize=(8,6))

for i in range(models.shape[0]):
    y = get_y(models.a[i], models.b[i])
    ax.plot(x, y, linestyle = 'dotted', color = 'grey')
y = get_y(ols['a'], ols['b'])

ax.plot(x, y, color='red', label='OLS')
ax.scatter(data.income, data.foodexp, alpha=.2)
ax.set_xlim((240, 3000))
ax.set_ylim((240, 2000))
legend = ax.legend()
ax.set_xlabel('Income', fontsize = 16)
ax.set_ylabel('Food Expenditure', fontsize = 16);<jupyter_output><empty_output><jupyter_text>Second plotThe 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.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.<jupyter_code>n = models.shape[0]
p1 = plt.plot(models.q, models.b, color='black', label='Quantile Reg.')
p2 = plt.plot(models.q, models.ub, color='black', linestyle='dotted')
p3 = plt.plot(models.q, models.lb, color='black', linestyle='dotted')
p4 = plt.plot(models.q, [ols['b']] * n, color='red', label='OLS')
P5 = plt.plot(models.q, [ols['lb']] * n, color='red', linestyle='dotted')
p6 = plt.plot(models.q, [ols['ub']] * n, color='red', linestyle='dotted')
plt.ylabel(r'$\beta_{income}$')
plt.xlabel('Quantiles of the conditional food expenditure distribution')
plt.legend()
plt.show();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15671222617487363, 'fraction_numerical': 0.10297697060475566, 'mean_word_length': 3.0408472012102874, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '349858', 'license': 'MIT', 'n_tokens_mistral': 2095, 'n_tokens_neox': 1743, 'n_words': 503, 'path': 'Quantile Regression.ipynb', 'repo_name': 'BeeChainFoundation/Brainstorming', 'size': '5341'}