Notebook 1: Empirical Relationships¶
Bryan Graham - University of California - Berkeley
Ec 240a: Econometrics, Fall 2015
The purpose of this iPython notebook, and the others which follow, is to concretely illustrate the application of some of the ideas introduced (more abstractly) in lecture. Python is a high-level general purpose computer programming language. It is freely available, and provides a convenient, and widely used, platform for scientific computing. The Anaconda distribution of Python 2.7 provides most of the core libraries we will use this semester and is easy to install (https://store.continuum.io/cshop/anaconda/).
There are many free resources available for learning Python. One, which may be of special interest to economists, can be found at http://quant-econ.net. A good introduction to data manipulation and organization in Python is provided by Wes McKinney in his book Python for Data Analysis.
Those components of course problem sets which involve computation should be completed in Python. Any exceptions to this requirement must be cleared with the course GSI. Course time will not be spent teaching Python. It is your responsibility to acquire the necessary coding skills independently outside of class. Your are, of course, free to ask for help from one another (as well as from me and the course GSI during scheduled office hours). I am mindful that students in Ec240a generally have a wide range of prior experiences with technical computing and programming. I ask that those of you with relatively strong backgrounds in these areas assist classmates who are less familiar with programming. I do not expect any of this to be easy, but I do hope that it is fun!
Python is one of a number of excellent technical computing environments. For most of my career I have used MATLAB, which I find very natural to code in, and for which many excellent high quality, commercial grade, libraries are available. MATLAB, along with many of its add-on libraries, is available to graduate students through the Econometrics Laboratory. I have also extensively used Mathematica at various times. Berkeley has a site license for this program. Both MATLAB and Mathematica are expensive, which may be a consideration. Another option is R, which is widely used in statistics. If speed is paramount, Julia may be of interest. I do not have much experience with these last two options, both of which, like Python, are free.
I strongly advise that you make a serious investment in at least one of these technical computing environments. You may be surprised at how helpful computing skills can be for various research tasks (not just econometrics).
Some of you may have prior experience with STATA. STATA is a commerical statistics and data analysis package. While STATA has a tremendous amount of functionality in certain areas, it is not a comprehensive technical computing environment (although the "MATA" part of STATA provides some functionality of this type). I frequently use STATA in my own work, and expect many of you will as well, but if your computation skills are confined to this program you will be limited as a researcher.
iPython Notebooks¶
This document is a static rendering (into html) of an iPython notebook. Mathematica, which some of you may be familiar with, also uses a "notebook" environment. Notebooks allow for the integration of text (including LaTex), graphics and code into a single document. They also provide a highly interactive environment for coding. For more complex programming tasks you might find it easier to write, debug and execute code in Spyder (which is included in the Anaconda distribution of Python). This provides a workspace environment similar to MATLAB or RStudio. However, for most of the assignments I will give in Ec240a the notebook approach will be both feasible and attractive.
# Direct Python to plot all figures inline (i.e., not in a separate window)
%matplotlib inline
# Load libraries
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import statsmodels.api as sm
# Division of two integers in Python 2.7 does not return a floating point result. The default is to round down
# to the nearest integer. The following piece of code changes the default.
from __future__ import division
Loading data¶
All the calculations in this notebook utilize the NLSY79_TeachingExtract.csv file which is available on the course webpage. This file includes variables extracted from multiple waves of the National Longitudinal Survey of Youth (1979). We will load this file into Python as a pandas dataframe. A good introduction to pandas is provided by the book Python for Data Analysis.
We will call our dataframe 'nlsy79' and use a hierarchical index structure. Rows are indexed by househould id (HHID_79) first, and then within households, by person id (PID_79).
# Directory where NLSY79_TeachingExtract.csv file is located
workdir = '/Users/bgraham/Dropbox/Teaching/Teaching_Datasets/'
# Read in NLSY79 Extract as a pandas dataframe
nlsy79 = pd.read_csv(workdir+'NLSY79_TeachingExtract.csv') # Reading .csv file as DataFrame
# Hierarchical index: household, then individual; keep indices as columns too
nlsy79.set_index(['HHID_79','PID_79'], drop=False)
nlsy79.rename(columns = {'AFQT_Adj':'AFQT'}, inplace=True) # Renaming AFQT
#Display the first few rows of the dataframe
nlsy79.head()
Our target sample will consist of all male, non-Black, non-Hispanic NLSY respondents. We will also exclude all units in the supplementary military and "poor White" subsamples. Follow-up of these two subsamples was discontinued early on. There are a total of 2,439 respondents in our target sample.
# Only retain non-black, non-hispanic, male NLSY79 respondents belonging to core sample
nlsy79 = nlsy79[(nlsy79.core_sample != 0) & (nlsy79.male != 0) & (nlsy79.black == 0) & (nlsy79.hispanic == 0)]
# Basic summary statistics of nlsy79 dataframe
nlsy79.describe()
In what follows we will explore the relationship between schooling and earnings as an adult. Our earnings variable will be the average of total annual earnings across the 1997, 1999, 2001 and 2003 calendar years. If earnings are missing for one of those years, the average will be over all non-missing years. NLSY79 respondents were born in the years 1957 to 1964 and hence were between the ages of 36 and 43 in 2000. Current earnings tend to be relatively close to permanent earnings at age 40, which motivates the selection of these calendar years. Earnings are measured in 2010 prices using the CPI-U price index.
In turns out that we are only able to measure earnings (and other key variables) for 1,969 of the respondents in our target sample of 2,439; that is for about 81% of target respondents. This missing data is due to a combination of item non-response and attrition. Appropriately addressing complications due to missing data is rather important in practice, but will be ignored here.
# Calculate average earnings across the 1997, 1999, 2001 and 2003 calendar years
# NOTE: This is an average over non-missing earnings values
nlsy79['Earnings'] = nlsy79[["real_earnings_1997", "real_earnings_1999", \
"real_earnings_2001", "real_earnings_2003"]].mean(axis=1)
# Only retain complete cases of year of birth, earnings, schooling, AFQT and parental education
nlsy79 = nlsy79[["PID_79","HHID_79","year_born","Earnings","HGC_Age28","AFQT"]]
nlsy79 = nlsy79.dropna()
# Summary statistics
nlsy79.describe()
Earnings is a continuously-valued variable (or practically so). For what follows it will be convenient to discretize earnings. The following snippet of code creates earnings buckets, in \$5,000 increments from \$0 to \$70,000, in \$10,000 increments from \$70,000 to \$100,000, and in \$50,000 increments from \$100,000 to \$250,000. All individuals with average earnings exceeding \$250,000 are lumped into a single bucket.
Schooling is discretely-valued, measured in whole years of completed schooling. However, in order to avoid complications due to sparse and/or empty cells, we will lump all individuals with less than 12 years of schooling (i.e., high school dropouts and less) into a single group.
After discretizing schooling and earnings we can compute the joint frequency distribution of schooling and earnings.
# Years-of-Schooling buckets
sch_bkts = [0, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]
# Annual earnings buckets
# NOTE: Buckets defined as monotonically decreasing so that the top row of the schooling/earnings
# cross-tabulation corresponds to the highest income bracket
earn_bkts = [2500000, 250000, 200000, 150000, 100000, 90000, 80000, 70000, 65000, 60000, 55000, 50000, \
45000, 40000, 35000, 30000, 25000, 20000, 15000, 10000, 5000, 0]
# Table of counts for each schooling-by-earnings cell with both row and column totals included
# NOTE: np.digitize() returns bin indices with upper bin edges excluded
earn_by_sch = pd.crosstab(np.digitize(nlsy79.Earnings, earn_bkts), np.digitize(nlsy79.HGC_Age28, sch_bkts), \
rownames=['Earnings'], colnames=['Schooling'], margins=True)
# Create informative labels for bins and attach them to the cross tabulation
sch_bin_labels = ['[0,12)', '12', '13', '14', '15', '16', '17', '18', '19', '20', 'Pr(Earn=y)']
earn_bin_labels = ['+250K','[200K,250K)','[150K,200K)','[100K,150K)','[90K,100K)','[80K,90K)','[70K,80K)', \
'[65K,70K)','[60K,65K)','[55K,60K)','[50K,55K)','[45K,50K)','[40K,45K)','[35K,40K)','[30K,35K)', \
'[25K,30K)','[20K,25K)','[15K,20K)','[10K,15K)','[5K,10K)','[0K,5K)','Pr(Sch=x)']
# Associate bin labels with cross tabulation
earn_by_sch.index = pd.Index(earn_bin_labels, name="Earnings")
earn_by_sch.columns = pd.Index(sch_bin_labels, name="Schooling")
# Print schooling-by-earnings joint frequency distribution (in percentages)
N = len(nlsy79.Earnings) # Total number of complete cases
print "Joint Frequency Distribution of Schooling and Earnings, NLSY79"
print ""
print np.round(100*earn_by_sch/N,2)
In inspecting the cross-tabulation we can see that most of the probability mass is concentrated along the diagonal running from low-schooling/low-earnings to high-schooling/high-earnings. However there is some mass is almost all cells.
Inspecting the marginal distribution of schooling (the bottom row of the table), we can see that there are large masses of probability in the schooling distribution at 12 and 16 years (corresponding to high school and college completion respectively).
The marginal distribution of earnings is reported in the right-most column of the table. Less than 3 percent of NLSY79 respondents in our subsample earned more than \$250,000 around 2000. About one half of respondents earned \$50,000 or less.
A somewhat more evocative way to report the joint frequency distribution is through a heatmap (with histograms depicting the two marginal distributions). The following snippet of code produces such a figure. In makes use of the seaborn library.
The code which follows also computes some basic probabilities and expectations. To do this we need to behave "as if" earnings is a discretely-valued random variable. For this purpose I assume that any individual in a given earnings bin earned the midpoint value of the bin. So I treat all individual in the \$0 to \$5,000 bin as earning \$2,500. I assume that all individuals in the \$250,000 and above bin earn \$350,000 per annum. This number was calibrated by making a Pareto Tail assumption.
# Convert cross tabulation into a pandas dataframe.
jfd = pd.DataFrame(100*earn_by_sch/N, index=earn_bin_labels, columns=sch_bin_labels)
# Calculate Pareto shape parameter for upper income bin
shp_num = np.log((np.array(earn_by_sch.loc[['+250K'],['Pr(Earn=y)']]) \
+ np.array(earn_by_sch.loc[['[200K,250K)'],['Pr(Earn=y)']])) \
/ np.array(earn_by_sch.loc[['[200K,250K)'],['Pr(Earn=y)']]))
shp_den = np.log(250000/200000)
alpha = shp_num/shp_den
# Attach numeric values to each earnings and schooling bin in order to calculate expectations
earn_lvls = np.array([round(float((alpha/(alpha-1))*250000) / 50000)*50000, 225000, 175000, 125000, 95000, 85000, 75000, \
67500, 62500, 57500, 52500, 47500, 42500, 37500, 32500, 27500, 22500, 17500, 12500, 7500, 2500])
sch_lvls = np.array([11, 12, 13, 14, 15, 16, 17, 18, 19, 20])
print ""
print "Calibrated Pareto scale parameter"
print float(alpha)
print ""
print "Calibrated average of upper earnings bin"
print round(float((alpha/(alpha-1))*250000) / 50000)*50000
# Drop row and columns corresponding to margins
jfd = jfd.drop('Pr(Sch=x)', 0) # Drops schooling margins
jfd = jfd.drop('Pr(Earn=y)', 1) # Drops earnings margins
# Get marginal distributions of earnings and schooling
mfd_earn = jfd.sum(axis=1)
mfd_sch = jfd.sum(axis=0)
# Calculate some probabilities and expectations
print ""
print "Probability of earning less than $50,000"
print jfd.loc[['[45K,50K)','[40K,45K)','[35K,40K)','[30K,35K)', \
'[25K,30K)','[20K,25K)','[15K,20K)','[10K,15K)','[5K,10K)','[0K,5K)']].sum(axis=1).sum(axis=0)
print ""
print "Probability of earning less than $50,000 and being a high school dropout"
print jfd.loc[['[45K,50K)','[40K,45K)','[35K,40K)','[30K,35K)', \
'[25K,30K)','[20K,25K)','[15K,20K)','[10K,15K)','[5K,10K)','[0K,5K)'],['[0,12)']].sum(axis=1).sum(axis=0)
print ""
print "Expected value of earnings"
print np.sum(mfd_earn/100 * earn_lvls)
print ""
print "Expected value of schooling"
print np.sum(mfd_sch/100 * sch_lvls)
# Form a plot with a heatmap of the joint frequency distribution at the center and histograms
# of the two marginal frequency distributions on each side.
import matplotlib.gridspec as gridspec
# Create figure; use gridspec to control relative size of three subplots
jfd_fig = plt.figure(figsize=(10, 10))
gs = gridspec.GridSpec(2, 2, width_ratios=[3, 1], height_ratios=[4,1])
# Joint frequency distribion
ax0 = plt.subplot(gs[0])
sns.set(style="white")
sns.heatmap(jfd, cmap="YlGnBu", robust=True, annot=True, fmt=".2f", cbar=False)
plt.title('Joint Frequency Distribution of Schooling and Earnings')
plt.xlabel('Years-of-Schooling')
plt.ylabel('Annual Earnings')
# Marginal frequency distribution of earnings
ax1 = plt.subplot(gs[1])
mfd_earn.plot(kind='barh',width=1, facecolor='blue', edgecolor='black').invert_yaxis()
plt.autoscale()
# Marginal frequency distribution of schooling
ax2 = plt.subplot(gs[2])
mfd_sch.plot(kind='bar',width=1, facecolor='blue', edgecolor='black')
plt.autoscale()
# Save result
plt.tight_layout()
plt.savefig(workdir+'Fig_jfd_earn_by_sch.png')
As discussed in lecture, we are often interested in the conditional distribution of one random variable given another. To get the conditional frequency distribution of earnings given schooling we divide each element of the joint frequency distribution by its column sum.
print "Conditional Frequency Distribution of Earnings Given Schooling, NLSY79"
print ""
cfd = (jfd / jfd.sum())*100
print np.round(cfd,2)
# Form a plot with a heatmap of the conditional frequency distribution at the center and histograms
# of the two marginal frequency distributions on each side.
# Create figure; use gridspec to control relative size of three subplots
cfd_fig = plt.figure(figsize=(10, 10))
gs = gridspec.GridSpec(2, 2, width_ratios=[3, 1], height_ratios=[4,1])
# Conditional frequency distribion
ax0 = plt.subplot(gs[0])
sns.set(style="white")
sns.heatmap(cfd, cmap="YlGnBu", robust=True, annot=True, fmt=".2f", cbar=False)
plt.title('Conditional Frequency Distribution of Earnings Given Schooling')
plt.xlabel('Years-of-Schooling')
plt.ylabel('Annual Earnings')
# Marginal frequency distribution of earnings
ax1 = plt.subplot(gs[1])
mfd_earn.plot(kind='barh',width=1, facecolor='blue', edgecolor='black').invert_yaxis()
plt.autoscale()
# Marginal frequency distribution of schooling
ax2 = plt.subplot(gs[2])
mfd_sch.plot(kind='bar',width=1, facecolor='blue', edgecolor='black')
plt.autoscale()
# Save result
plt.tight_layout()
plt.savefig(workdir+'Fig_cfd_earn_given_sch.png')
Conditional mean and median regression¶
The conditional distribution of earnings given schooling is easier to interpret than the corresponding joint one. For example, we observe that while 25 percent of respondents with 20 years of schooling (essentially corresponding to a Ph.D) earned more that $250,000, only 0.5 percent of high school graduates did. The heatmap also clearly shows that the conditional distribution of earnings "shifts to the right" as schooling increases.
To get a single "prototypical" measure of earnings for each schooling level, we can use the pandas groupby function to compute mean and median earnings conditional on schooling.
# Assign NLSY respondents to a discrete schooling group
nlsy79['Discrete_Schooling'] = np.digitize(nlsy79.HGC_Age28, sch_bkts)
# Group nlsy79 dataframe by schooling group
grouped_earn = nlsy79['Earnings'].groupby(nlsy79['Discrete_Schooling'])
# Create a dataframe with each row corresponding to a schooling group
regframe = pd.DataFrame({'Schooling': sch_bin_labels[0:10]},index=[1,2,3,4,5,6,7,8,9,10])
# Calculate mean, median and the 0.25 and 0.75 quantiles of earnings for each schooling group
regframe['Mean'] = grouped_earn.mean()
regframe['Std'] = grouped_earn.std()
regframe['Mean_plus_1Std'] = regframe['Mean'] + regframe['Std']
regframe['Mean_minus_1Std'] = regframe['Mean'] - regframe['Std']
regframe['Median'] = grouped_earn.median()
regframe['0.25 Quantile'] = grouped_earn.quantile(q=0.25)
regframe['0.75 Quantile'] = grouped_earn.quantile(q=0.75)
print "Mean, standard deviation, median and the 0.25 and 0.75 quantiles of earnings for each schooling group"
print ""
print regframe
# Plot mean earnings +/- 1 std by schooling group
reg_fig = plt.figure(figsize=(10, 10))
ax1 = reg_fig.add_subplot(1,2,1)
sns.pointplot(x="Schooling", y="Mean", data=regframe, color="black", marker="o")
sns.pointplot(x="Schooling", y="Mean_plus_1Std", data=regframe, color="gray", markers=["x"], linestyles=["--"])
sns.pointplot(x="Schooling", y="Mean_minus_1Std", data=regframe, color="gray", markers=["x"], linestyles=["--"])
plt.title('Mean & +/- 1 Std. Dev.')
plt.xlabel('Years-of-Schooling')
plt.ylabel('Annual Earnings')
plt.ylim([0,250000])
# Plot median earnings and interquartile range of earnings by schooling group
ax2 = reg_fig.add_subplot(1,2,2)
sns.pointplot(x="Schooling", y="0.25 Quantile", data=regframe, color="gray", \
markers=["x"], linestyles=["--"])
sns.pointplot(x="Schooling", y="Median", data=regframe, color="black")
sns.pointplot(x="Schooling", y="0.75 Quantile", data=regframe, color="gray", \
markers=["x"], linestyles=["--"])
plt.title('Median & Interquartile Range')
plt.xlabel('Years-of-Schooling')
plt.ylabel('Annual Earnings')
plt.ylim([0,250000])
# Save plot
plt.tight_layout()
plt.savefig(workdir+'Fig_cef_earn_given_sch.png')
We will use this same subsample of NLSY79 respondents in future notebooks. We therefore save select variables in our nlsy79 dataframe to a comma delimited file for future use.
# Extract portion of dataframe to retain
nlsy79 = nlsy79[["PID_79","HHID_79","year_born","Earnings","HGC_Age28","AFQT"]]
# Export estimation subsample to comma delimited file for later use
nlsy79.to_csv(path_or_buf=workdir+'nlsy79.csv')
# This imports an attractive notebook style from Github
from IPython.core.display import HTML
import urllib2
HTML(urllib2.urlopen('http://bit.ly/1Bf5Hft').read())