sales-hedonic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial on Hedonic Regression\n",
    "\n",
    "This material uses Python to demonstrate some aspects of hedonic regression, but the objective here is not to learn to program, but to understand the hedonic regression methodology.  We begin with an example in which we generate some synthetic data using a set of coefficients and a mathematical model, and learn those coefficients using a statistical method called multiple regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Startup steps\n",
    "import pandas as pd, numpy as np, statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt, matplotlib.cm as cm, matplotlib.font_manager as fm\n",
    "import matplotlib.mlab as mlab\n",
    "from scipy.stats import pearsonr, ttest_rel\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using a model we define:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nsample = 100\n",
    "X = np.linspace(0, 10, 100)\n",
    "beta = np.array([0, 2])\n",
    "e = np.random.normal(size=nsample)\n",
    "X = sm.add_constant(X)\n",
    "y = np.dot(X, beta) + e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the data and the model.  Note that the intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(1, figsize=(10,8), )\n",
    "plt.plot([0, 10], [0, 20])\n",
    "plt.scatter(x, y, marker=0, s=10, c='g')\n",
    "plt.axis([0, 10, 0, 20])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = sm.OLS(y, X)\n",
    "results = model.fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('Parameters: ', results.params)\n",
    "print('R2: ', results.rsquared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nsample = 100\n",
    "x = np.linspace(0, 10, 100)\n",
    "X = np.column_stack((x, x**2))\n",
    "beta = np.array([1, 2, .5])\n",
    "e = np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = sm.add_constant(X)\n",
    "y = np.dot(X, beta) + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(x,y, s=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = sm.OLS(y, X)\n",
    "results = model.fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sf = pd.read_csv('data/redfin_2017-03-05-17-45-34-san-francisco-county-1-month.csv')\n",
    "sf.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sf1 = sf.rename(index=str, columns={'SALE TYPE': 'saletype',\n",
    "    'SOLD DATE': 'solddate', 'PROPERTY TYPE': 'proptype', 'ADDRESS': 'address',\n",
    "    'CITY': 'city', 'STATE': 'state', 'ZIP': 'zip', 'PRICE': 'price', 'BEDS': 'beds',\n",
    "    'BATHS': 'baths', 'LOCATION': 'location', 'SQUARE FEET': 'sqft', 'LOT SIZE': 'lotsize',\n",
    "    'YEAR BUILT': 'yrbuilt', 'DAYS ON MARKET': 'daysonmkt', '$/SQUARE FEET': 'pricesqft',\n",
    "    'LATITUDE': 'latitude', 'LONGITUDE': 'longitude', 'HOA/MONTH': 'hoamonth',\n",
    "    'URL (SEE http://www.redfin.com/buy-a-home/comparative-market-analysis FOR INFO ON PRICING)': 'url',\n",
    "    'STATUS': 'status', 'NEXT OPEN HOUSE START TIME': 'nextopenstart', 'NEXT OPEN HOUSE END TIME': 'nextopenend',\n",
    "    'SOURCE': 'source', 'MLS#': 'mls', 'FAVORITE': 'favorite', 'INTERESTED': 'interested'\n",
    "    })\n",
    "\n",
    "sf1.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sf1.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(1, figsize=(10,8), )\n",
    "plt.scatter(sf1['sqft'], sf1['price'], marker=0, s=10, c='g')\n",
    "#plt.axis([12, 16, 12, 16])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import numpy as np\n",
    "from patsy import dmatrices\n",
    "y, X = dmatrices('np.log(price) ~ np.log(sqft) + C(baths)', \n",
    "                 data=sf1, return_type='dataframe')\n",
    "mod = sm.OLS(y, X)\n",
    "res = mod.fit()\n",
    "residuals = res.resid\n",
    "predicted = res.fittedvalues\n",
    "observed = y\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(residuals, bins=25, normed=True, alpha=.5)\n",
    "mu = residuals.mean()\n",
    "variance = residuals.var()\n",
    "sigma = residuals.std()\n",
    "x = np.linspace(-3, 3, 100)\n",
    "plt.plot(x,mlab.normpdf(x, mu, sigma));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(1, figsize=(10,8), )\n",
    "plt.plot([12, 16], [0, 0], c='b')\n",
    "plt.scatter(predicted, residuals, marker=0, s=10, c='g');\n",
    "plt.axis([12, 16, -0.8, 0.8])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(1, figsize=(10,8), )\n",
    "plt.plot([12, 16], [12, 16])\n",
    "plt.scatter(observed, predicted, marker=0, s=10, c='g')\n",
    "plt.axis([12, 16, 12, 16])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vars = pd.read_csv('data/ba_block_variables.csv', dtype={'block_id':object})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vars.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "smallvars = vars[['block_id','block_groups_sum_persons','block_groups_sum_acres','block_groups_total_jobs',\n",
    "                 'block_groups_median_children', 'block_groups_median_persons' , 'block_groups_median_income',\n",
    "                 'block_groups_prop_tenure_2', 'nodes_low_income_hh_1500m', 'nodes_high_income_hh_1500m',\n",
    "                 'nodes_jobs_5000m', 'nodes_jobs_30km', 'nodes_population_400m', 'nodes_population_800m',\n",
    "                 'nodes_jobs_800m', 'prop_race_of_head_1', 'prop_race_of_head_2', 'pumas_density_households',\n",
    "                 'nodes_jobs_3000m_agg1', 'nodes_jobs_3000m_agg2', 'pumas_density_jobs', 'nodes_jobs_40km',\n",
    "                 'nodes_jobs_3000m_agg3', 'nodes_jobs_3000m_agg4', 'nodes_jobs_3000m_agg5',\n",
    "                 'block_groups_prop_persons_1', 'block_groups_prop_persons_2', 'block_groups_median_age_of_head',\n",
    "                 'puma10_id_is_0607501', 'puma10_id_is_0607502', 'puma10_id_is_0607503', 'puma10_id_is_0607504']]\n",
    "sv = smallvars.rename(index=str, columns={'block_groups_sum_persons': 'bgpop',\n",
    "                                     'block_groups_sum_acres': 'bgacres',\n",
    "                                    'block_groups_total_jobs': 'bgjobs',\n",
    "                                    'block_groups_median_children': 'bgmedkids',\n",
    "                                    'block_groups_median_persons': 'bgmedhhs',\n",
    "                                    'block_groups_median_income': 'bgmedinc',\n",
    "                                    'block_groups_prop_tenure_2': 'proprent',\n",
    "                                    'nodes_low_income_hh_1500m': 'lowinc1500m',\n",
    "                                    'nodes_high_income_hh_1500m': 'highinc1500m',\n",
    "                                    'nodes_jobs_5000m': 'lnjobs5000m',\n",
    "                                    'nodes_jobs_30km': 'lnjobs30km',\n",
    "                                    'nodes_jobs_40km': 'lnjobs40km',\n",
    "                                    'nodes_population_400m': 'lnpop400m',\n",
    "                                    'nodes_population_800m': 'lnpop800m',\n",
    "                                    'nodes_jobs_800m': 'lnjobs800m',\n",
    "                                    'prop_race_of_head_1': 'propwhite',\n",
    "                                    'prop_race_of_head_2': 'propblack',\n",
    "                                    'pumas_density_households': 'pumahhden',\n",
    "                                    'pumas_density_jobs': 'pumajobden',\n",
    "                                    'nodes_jobs_3000m_agg1': 'lnbasic3000m',\n",
    "                                    'nodes_jobs_3000m_agg2': 'lntcpuw3000m',\n",
    "                                    'nodes_jobs_3000m_agg3': 'lnret3000m',\n",
    "                                    'nodes_jobs_3000m_agg4': 'lnfire3000m',\n",
    "                                    'nodes_jobs_3000m_agg5': 'lnserv3000m',\n",
    "                                    'block_groups_prop_persons_1': 'prop1per',\n",
    "                                    'block_groups_prop_persons_2': 'prop2per',\n",
    "                                    'block_groups_median_age_of_head': 'bgmedagehd',\n",
    "                                    'puma10_id_is_0607501': 'puma1',\n",
    "                                    'puma10_id_is_0607502': 'puma2',\n",
    "                                    'puma10_id_is_0607503': 'puma3',\n",
    "                                    'puma10_id_is_0607504': 'puma4'\n",
    "                                         })\n",
    "sv.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recoded detailed race code\n",
    " 1 .White alone\n",
    " 2 .Black or African American alone\n",
    " 3 .American Indian alone\n",
    " 4 .Alaska Native alone\n",
    " 5 .American Indian and Alaska Native tribes specified; or American\n",
    " .Indian or Alaska native, not specified and no other races\n",
    " 6 .Asian alone\n",
    " 7 .Native Hawaiian and Other Pacific Islander alone\n",
    " 8 .Some other race alone\n",
    " 9 .Two or more major race groups "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open('pg_engine_redfin.txt') as f:\n",
    "    pg_engine = f.readlines()\n",
    "from sqlalchemy import create_engine\n",
    "engine = create_engine(pg_engine[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "import pandas as pd\n",
    "df = pd.read_sql_query('select * from \"rental_listings\"',con=engine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(df.dtypes)\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# convert the date column to yyyy-mm-dd date format\n",
    "df['date'] = pd.to_datetime(df['date'], format='%Y-%m-%d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# load the 2014 census data set of MSAs\n",
    "census = pd.read_csv('data/census_pop_income.csv', encoding='ISO-8859-1')\n",
    "census['median_income'] = census['2014_median_income'].str.replace(',','').astype(int)\n",
    "census['population'] = census['2014_pop_est'].str.replace(',','').astype(int)\n",
    "census = census.drop(labels='notes', axis=1, inplace=False)\n",
    "census = census.set_index('region')\n",
    "census.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# these are the 15 most populous metros by population, defined by census bureau 2014 estimates\n",
    "most_populous_regions = census['2014_pop_est'].sort_values(ascending=False, inplace=False)\n",
    "print(most_populous_regions.head(15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = df.merge(census, left_on='region', right_index=True)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create an HDF5 file if desired, in the data directory\n",
    "#df.to_hdf('data/rents.h5','rents',append=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load from HDF5 if desired\n",
    "#df = pd.HDFStore('data/rents.h5')\n",
    "#df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "upper_percentile = 0.998\n",
    "lower_percentile = 0.002\n",
    "\n",
    "# how many rows would be within the upper and lower percentiles?\n",
    "upper = int(len(df) * upper_percentile)\n",
    "lower = int(len(df) * lower_percentile)\n",
    "\n",
    "# get the rent/sqft values at the upper and lower percentiles\n",
    "rent_sqft_sorted = df['rent_sqft'].sort_values(ascending=True, inplace=False)\n",
    "upper_rent_sqft = rent_sqft_sorted.iloc[upper]\n",
    "lower_rent_sqft = rent_sqft_sorted.iloc[lower]\n",
    "\n",
    "# get the rent values at the upper and lower percentiles\n",
    "rent_sorted = df['rent'].sort_values(ascending=True, inplace=False)\n",
    "upper_rent = rent_sorted.iloc[upper]\n",
    "lower_rent = rent_sorted.iloc[lower]\n",
    "\n",
    "# get the sqft values at the upper and lower percentiles\n",
    "sqft_sorted = df['sqft'].sort_values(ascending=True, inplace=False)\n",
    "upper_sqft = sqft_sorted.iloc[upper]\n",
    "lower_sqft = sqft_sorted.iloc[lower]\n",
    "\n",
    "print('valid rent_sqft range:', [lower_rent_sqft, upper_rent_sqft])\n",
    "print('valid rent range:', [lower_rent, upper_rent])\n",
    "print('valid sqft range:', [lower_sqft, upper_sqft])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create a boolean vector mask to filter out any rows with rent_sqft outside of the reasonable values\n",
    "rent_sqft_mask = (df['rent_sqft'] > lower_rent_sqft) & (df['rent_sqft'] < upper_rent_sqft)\n",
    "\n",
    "# create boolean vector masks to filter out any rows with rent or sqft outside of the reasonable values\n",
    "rent_mask = (df['rent'] > lower_rent) & (df['rent'] < upper_rent)\n",
    "sqft_mask = (df['sqft'] > lower_sqft) & (df['sqft'] < upper_sqft)\n",
    "\n",
    "# filter the thorough listings according to these masks\n",
    "filtered_listings = pd.DataFrame(df[rent_sqft_mask & rent_mask & sqft_mask])\n",
    "len(filtered_listings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filtered_listings.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#sfbay = filtered_listings[filtered_listings['region']=='sfbay']\n",
    "sfbay = df[df['region']=='sfbay']\n",
    "sfbay.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#sfbay['rent_sqft'].quantile(.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#sfbay['sqft'].quantile(.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create a boolean vector mask to filter out any rows with rent_sqft and sqft in Bay Area under 1 percentile\n",
    "#sfbay_rent_sqft_mask = (sfbay['rent_sqft'] > sfbay['rent_sqft'].quantile(.01) )\n",
    "\n",
    "# create boolean vector masks to filter out any rows with rent or sqft outside of the reasonable values\n",
    "#sfbay_sqft_mask = (sfbay['sqft'] > sfbay['sqft'].quantile(.01) )\n",
    "\n",
    "# filter the thorough listings according to these masks\n",
    "#sfbay_filtered = pd.DataFrame(sfbay[sfbay_rent_sqft_mask & sfbay_sqft_mask])\n",
    "#len(sfbay_filtered)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "upper_percentile = 0.998\n",
    "lower_percentile = 0.002\n",
    "\n",
    "# how many rows would be within the upper and lower percentiles?\n",
    "upper = int(len(sfbay) * upper_percentile)\n",
    "lower = int(len(sfbay) * lower_percentile)\n",
    "\n",
    "# get the rent/sqft values at the upper and lower percentiles\n",
    "rent_sqft_sorted = sfbay['rent_sqft'].sort_values(ascending=True, inplace=False)\n",
    "upper_rent_sqft = rent_sqft_sorted.iloc[upper]\n",
    "lower_rent_sqft = rent_sqft_sorted.iloc[lower]\n",
    "\n",
    "# get the rent values at the upper and lower percentiles\n",
    "rent_sorted = sfbay['rent'].sort_values(ascending=True, inplace=False)\n",
    "upper_rent = rent_sorted.iloc[upper]\n",
    "lower_rent = rent_sorted.iloc[lower]\n",
    "\n",
    "# get the sqft values at the upper and lower percentiles\n",
    "sqft_sorted = sfbay['sqft'].sort_values(ascending=True, inplace=False)\n",
    "upper_sqft = sqft_sorted.iloc[upper]\n",
    "lower_sqft = sqft_sorted.iloc[lower]\n",
    "\n",
    "print('valid rent_sqft range:', [lower_rent_sqft, upper_rent_sqft])\n",
    "print('valid rent range:', [lower_rent, upper_rent])\n",
    "print('valid sqft range:', [lower_sqft, upper_sqft])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create a boolean vector mask to filter out any rows with rent_sqft outside of the reasonable values\n",
    "rent_sqft_mask = (sfbay['rent_sqft'] > lower_rent_sqft) & (sfbay['rent_sqft'] < upper_rent_sqft)\n",
    "\n",
    "# create boolean vector masks to filter out any rows with rent or sqft outside of the reasonable values\n",
    "rent_mask = (sfbay['rent'] > lower_rent) & (sfbay['rent'] < upper_rent)\n",
    "sqft_mask = (sfbay['sqft'] > lower_sqft) & (sfbay['sqft'] < upper_sqft)\n",
    "\n",
    "# filter the thorough listings according to these masks\n",
    "sfbay_filtered = pd.DataFrame(sfbay[rent_sqft_mask & rent_mask & sqft_mask])\n",
    "len(sfbay_filtered)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sfbay_filtered.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sf = sfbay_filtered.merge(sv, left_on='fips_block', right_on='block_id')\n",
    "sf['northsf'] = sf['puma1']+sf['puma2']+sf['puma3']+sf['puma4']\n",
    "sf['sqft1'] = sf['sqft']*sf['sqft']<1500\n",
    "sf['sqft2'] = sf['sqft']*sf['sqft']>1499\n",
    "sf['pct1per'] = sf['prop1per']*100\n",
    "sf['pct2per'] = sf['prop2per']*100\n",
    "sf['pctrent'] = sf['proprent']*100\n",
    "sf['pctblack'] = sf['propblack']*100\n",
    "sf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sf.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sf['bgpopden'] = sf['bgpop']/sf['bgacres']\n",
    "sf['bgjobden'] = sf['bgjobs']/sf['bgacres']\n",
    "sf['highlowinc1500m'] = sf['highinc1500m']/sf['lowinc1500m']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sf.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use either filtered_listings for the national dataset or sfbay_filtered for the Bay Area subset\n",
    "\n",
    "import statsmodels.api as sm\n",
    "import numpy as np\n",
    "from patsy import dmatrices\n",
    "y, X = dmatrices('np.log(rent) ~ np.log(sqft) + bedrooms + bathrooms \\\n",
    "                 ', \n",
    "                 data=sf, return_type='dataframe')\n",
    "mod = sm.OLS(y, X)\n",
    "res = mod.fit()\n",
    "residuals = res.resid\n",
    "predicted = res.fittedvalues\n",
    "observed = y\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use either filtered_listings for the national dataset or sfbay_filtered for the Bay Area subset\n",
    "\n",
    "import statsmodels.api as sm\n",
    "import numpy as np\n",
    "from patsy import dmatrices\n",
    "y, X = dmatrices('np.log(rent) ~ np.log(sqft) + bedrooms + bathrooms  +\\\n",
    "                 np.log(bgmedinc) + lnjobs5000m + lnjobs40km + pctblack \\\n",
    "                  + pumahhden + pctrent + pct1per + pct2per + bgmedagehd  \\\n",
    "                 ', \n",
    "                 data=sf, return_type='dataframe')\n",
    "mod = sm.OLS(y, X)\n",
    "res = mod.fit()\n",
    "residuals = res.resid\n",
    "predicted = res.fittedvalues\n",
    "observed = y\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(residuals, bins=25, normed=True, alpha=.5)\n",
    "mu = residuals.mean()\n",
    "variance = residuals.var()\n",
    "sigma = residuals.std()\n",
    "x = np.linspace(-3, 3, 100)\n",
    "plt.plot(x,mlab.normpdf(x, mu, sigma));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(1, figsize=(10,8), )\n",
    "plt.plot([7, 9], [0, 0], c='b')\n",
    "plt.scatter(predicted, residuals, marker=0, s=2, c='g');\n",
    "plt.axis([7.25, 9, -1.5, 1.5])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(1, figsize=(10,8), )\n",
    "plt.plot([6, 9.5], [6, 9.5])\n",
    "plt.scatter(observed, predicted, marker=0, s=2, c='g')\n",
    "plt.axis([6.5, 9.5, 6.5, 9.5])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(residuals.mean())\n",
    "print(residuals.std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# If we want to use WLS we need a useful set of weights. The default produces the same results as OLS \n",
    "mod = sm.WLS(y, X, weights=1.)\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Warning, this is a very intensive process and will take a while\n",
    "%%time\n",
    "from pymc3 import Model, NUTS, sample\n",
    "from pymc3.glm import glm\n",
    "\n",
    "with Model() as model_glm:\n",
    "    glm('np.log(rent) ~ np.log(sqft) + bedrooms + bathrooms', sfbay_filtered)\n",
    "    trace = sample(5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from pymc3 import traceplot\n",
    "%matplotlib inline\n",
    "traceplot(trace);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import optimize\n",
    "from pymc3 import find_MAP\n",
    "with model_glm:\n",
    "\n",
    "    # obtain starting values via MAP\n",
    "    start = find_MAP(fmin=optimize.fmin_powell)\n",
    "\n",
    "    # draw 2000 posterior samples\n",
    "    trace = sample(5000, start=start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "traceplot(trace);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import theano\n",
    "import pymc3 as pm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(7, 7))\n",
    "ax = fig.add_subplot(111, xlabel='x', ylabel='y', title='Generated data and underlying model')\n",
    "ax.plot(np.log(sfbay['sqft']), np.log(sfbay['rent']), 'o', markersize=.5, color='blue', label='sampled data')\n",
    "#ax.plot(x, true_regression_line, label='true regression line', lw=2.)\n",
    "#pm.glm.plot_posterior_predictive(trace, samples=100,\n",
    "#                                 label='posterior predictive regression lines')\n",
    "plt.legend(loc=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}